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A3STRACT 

The  basic  principles  of  the  Galericin  finite  element 
method  are  discussed  and  applied  to  two  different 
formulations;  one  using  different  basis  functions  and  the 
other  using  the  vorticity-divergence  form  of  the  shallow 
water  equations.  Each  formulation  is  compared  to  the 
primitive  form  of  the  equations  developed  by  lelley  (1976). 
The  testing  involves  a  comparison  of  three  finite  element 
prediction  models  using  variable  size  elements.  Equilateral 
elements  significantly  improve  the  solution  and  are  used  in 
most  of  the  comparisons.  The  formulation  using  different 
basis  functions  produces  poorer  results  than  the  primitive 
formulation .  The  vorticity-divergence  formulation  produces 
superior  results  while  executing  faster  than  the  primitive 
model.  Fcwever,  it  does  require  more  storage  and  the 
relaxation  parameters  are  sensitive  to  the  domain  geometry. 
The  computer  implementation  for  the  vorticity-divergence 
model  is  discussed  and  the  source  listing  is  included. 
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I.  INTROEUCTION 


Shuman  fl978)  claims  that  progress  in  numerical  modeling 
of  the  general  circulation  has  been  to  some  degree  dictated 
In  the  past  by  the  rate  of  development  In  the  field  of 
computer  technology.  However,  the  limited  ability  to 
oarameterlze  the  effects  of  small-scale  processes  In  terms 
of  large  scale  motions  has  been  an  equally  Important 
limiting  factor.  Essentially,  the  major  problem  of  numerical 
modeling  of  the  general  circulation  Is  simply  that  of 
producing  a  very  long  range  numerical  weather  forecast. 

Certainly  the  equations  used  In  the  models  must  be  more 
sophisticated  to  include  those  physical  processes  which-  are 
unimportant  for  a  short  range  forecast,  but  may  become 
crucial  as  the  length  of  the  forecast  Is  extended.  Another 
area  where  concentrated  efforts  have  Improved  the  forecast 
Involves  the  computational  techniques  employed  to 
approximate  and  solve  the  governing  equations  of  the  models. 

The  motivation  behind  this  thesis  is  to  Investigate  the 
application  of  a  relatively  new  computational  technique  to 
the  field  of  numerical  weather  prediction.  The  finite 
element  method,  long  established  in  engineering,  has  been 
seriously  considered  only  during  the  past  decade  In 
meteorology.  This  method  has  great  potential  for  application 
In  atmospheric  prediction  models. 
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A.  BACKGROUND 


The  most  common  numerical  integration  procedure  for 
weather  prediction  has  been  the  finite  difference  method  in 
which  the  derivatives  in  the  differential  equations  of 
motion  are  replaced  by  finite  difference  approximations  at  a 
discrete  set  of  points  in  space  and  time.  The  resulting  set 
of  equations,  with  appropriate  restrictions,  can  then  he 
solved  by  algebraic  methods.  Until  recently,  the  finite 
difference  method  has  been  the  workhorse  in  atmospheric 
prediction  models,  from  their  first  computer  implementation 
to  the  present. 

With  the  introduction  of  each  new  generation  of 
computers.  the  gap  between  numerical  forecasts  and 
atmospheric  observations  has  decreased.  The  rate  at  which 
this  gap  decreased  has  slowed  down  and  appears  to  be 
leveling  off.  This  would  Indicate  that  computer  technology 
may  not  be  the  primary  obstruction  to  better  numerical 
forecasts.  In  fact,  bigger  and  faster  computers  alone  have 
demonstrated  their  inability  to  significantly  Improve  the 
numerical  forecast. 

For  example,  a  major  limiting  factor  of  finite 
difference  approximations  is  the  truncation  error.  The 
National  Weather  Service  7  Layer  Primitive  Equation  model 
(7LPE  Model),  operational  from  1966  to  1980,  had  truncation 
errors  which  increased  at  a  rate  proportional  to  the  square 
o*  the  grid  spacing.  That  is,  the  smaller  the  grid  interval. 


the  smaller  the  truncation  error.  To  increase  its  accuracy 
would  require  increasing  the  grid  matrix  density.  This  would 
require  increased  computer  storage  and  computational  time. 
State  of  the  art  computers  are  capable  of  providing  these 
additional  resources. 

The  problem  now  goes  beyond  numerical  techniques  and 
computer  technology.  Operationally,  the  National  Weather 
service  is  not  capable  (due  to  monetary  restrictions)  of 
providing  a  denser  concentration  of  atmospheric 
observations.  Therefore,  with  the  present  density  of  Initial 
data  (observations)  and  objective  analysis  techniques 
(getting  the  data  for  grid  points  by  interpolating  from 
observed  data  sources),  reducing  the  grid  spacing  further  on 
the  71PS  Mode]  does  not  significantly  increase  the  accuracy 
cf  the  solution. 

This  additional  computer  capability  can  not  be  utilized 
using  finite  difference  methods.  Therefore,  new  numerical 
integration  techniques  must  be  investigated,,  such  that  given 
the  same  density  of  observed  lata,  superior  solutions  are 
produced. 

Two  alternative  techniques,  the  spectral  method  and  the 
finite  element  method,  have  started  to  gain  attention.  Both 
the  spectral  and  finite  element  methods  require  more 
computational  time  per  forecast  time  step  than  does  the 
finite  difference  method.  For  example,  the  finite  element 
method  requires  an  equation  solver  to  invert  a  larger  matrix 


at  each  tire  step  for  each  variable.  In  this  sease,  these 
methods  were  held  back  by  computer  technology,  but  recent 
advances  in  computer  technology  (i.e.  larger  and  faster 
storage  devices,  multi-processors,  etc)  have  made  these 
alternative  numerical  techniques  competitive. 

For  long  range  weather  predictions,  the  spectral  method 
applied  over  the  globe  or  hemisphere  is  a  natural  method, 
due  to  the  existence  of  efficient  transforms  for  the 
nonlinear  terms  on  spherical  geometry.  It  also  eliminates 
the  truncation  error  for  the  horizontal  space  derivatives 
and  the  nonlinear  instability  (aliasing).  For  these  reasons, 
global  spectral  models  have  been  developed  and  implemented 
on  an  operational  level,  replacing  the  global  finite 
difference  models. 

However,  because  the  spectral  harmonics  are  globally 
rather  than  locally  defined,  it  is  thought  that  for  problems 
of  more  detailed  limited  area  forecasting,  the  finite 
element  method  is  more  suitable.  Pioneering  voxk  to  adapt 
finite  element  methods  to  meteorological  applications  has 
been  done  by  Cullen  (1973, 1974  and  1S79),  Staniforth  and 
Mitchell  C1977),  Hinsman  (1975)  and  Kelley  (1976).  The  most 
recent  finite  element  meteorological  model  at  the  Naval 
Postgraduate  School  was  written  by  Kelley  (1976)  with  the 
collaboration  of  Ir.  R.T.  Williams.  It  is  this  study  that 
will  serve  as  a  basis  for  this  thesis.  The  model  written  by 
Kelley  will  be  altered  and  used  for  comparative  testing  with 
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improved  finite  element  forms  implemented  by  this  author. 
Some  of  the  techniques  and  coles  developed  by  Kelley  are 
alsc  employed  in  this  thesis.  Older  (1981)  developed  a 
technique  to  smoothly  vary  the  grid  geometry  in  the  domain. 
This  technique  is  also  implemented  both  on  Kelley's  model 
and  with  the  new  formulation  to  give  greater  versatility 
when  testing  the  model  performance. 

3.  C3JECTI7ES 

The  objectives  for  this  thesis  can  be  divided  into  two 
categories:  1)  meteorology,  2)  computer  science.  First,  the 
meteorological  objectives  of  developing  improved  finite 
element  forms  for  shallow  water  equations  are  as  follows: 

1)  -  Older  (1981)  after  collaboration  with  Dr. 
M.J.P.  Cullen,  showed  how  equilaterally  shaped  elements 
produced  significantly  better  results  than  did  other 
triangular  elements.  Kelley  (1976)  used  right  triangular 
elements  in  the  implementation  of  a  two-dimensional  finite 
element  model  using  the  primitive  form  of  the  shallow  water 
equations.  A  considerable  amount  of  small-scale  noise  was 
observed  in  the  solution.  Hereafter,  this  model,  which  was 
developed  by  Kelley  (1976),  will  be  referred  to  as  the 
primitive  model.  This  first  objective  involves 
re-implementing  the  primitive  model  using  equilaterally 
shaped  elements  and  comparing  the  results  to  those  in 
Kelley's  thesis. 
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2)  -  Williams  and  Zienkiewicz  (19eD  presented  new 
finite  element  techniques  for  formulations  for  the  shallow 
water  eolations,  which  use  differently  shaped  functions  to 
approximate  the  different  dependent  variables,  which  in 
effect  stagger  the  variab^^s.  Schoenstadt  (1980) 
demonstrated  the  advantage  of  spatial  staggering  of 
dependent  variables  in  finite  difference  models.  The 
application  of  this  technique  to  finite  element  models  is  a 
natural  extension,  and  excellent  results  were  obtained  by 
Williams  and  Zienklewicz  (1981)  from  application  of  these 
new  formulations  on  linearized  one  dimensional  cases.  The 
objective  here  is  to  implement  the  new  forms  on  the 
primitive  model  and  again  do  quanitative  comparisons  of  the 
results. 

7)  -  The  major  emphasis  in  this  study  deals  with  the 
implementation  and  comparison  of  the  vortlcity  divergence 
form  of  the  shallow  water  equations,  which  is  described  in 
detail  in  Chapter  HI.  This  formulation  has  the  following 
advantages.  First,  the  geostrophlc  adjustment  process  is 
treated  better  than  la  the  primitive  form  of  the  equations. 
Secondly,  the  velocity  and  height  fields  are  evaluated  at 
the  same  grid  point,  where  the  best  primitive  form  requires 
staggering  these  dependent  variables.  And  thirdly,  a  larger 
time  step  is  allowed  due  to  the  semi  implicit  form  of  the 
equation.  Again  comparisons  between  the  results  from  the 
vortlcity  divergence  and  primitive  model  are  presented. 
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The  computer  science  aspect  of  this  thesis  was  primarily 
devoted  to  the  Implementation  of  the  different  models  and 
the  style  and  architecture  of  the  program.  Finite  element 
methods  require  more  computational  time  than  do  finite 
difference  methods,  not  only  in  the  solution  of  the 
equations,  hut  also  in  the  amount  of  computation  required  to 
evaluate  each  term  In  the  equations. 

The  Implementations  of  these  two  dimensional  models, 
although  complex  when  viewed  from  the  surface,  have  a  lot  of 
generality  and  redundancy  In  the  operations  required. 
Versatile  modules  can  be  written  to  ease  the  Implementation 
and  facilitate  changes.  The  objective  here  Is  to  efficiently 
Implement  these  new  forms  and  demonstrate  the  utility  of 
these  versatile  modules  for  future  Implementations. 

C.  THESIS  STRUCTURE 

This  thesis  presents  the  results  obtained  from  tests  of 
the  various  finite  element  .  formulations.  The  results  are 
compared  to  those  from  the  primitive  model  written  by  Selley 
(19?6).  Accompanying  the  results  Is  a  detailed  discussion  of 
the  reformulation  and  Implementation  process. 

Chapter  II  of  the  thesis  presents  a  tutorial  of  the 
finite  element  method  and  the  area  coordinates  system  used 
In  the  evaluation  of  the  element  integration.  The  Galericln 
finite  element  method  used  in  all  the  models  Is  developed 
and  applied  to  the  advection  equation  In  one  dimension. 
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Chapter  III  presents  the  detailed  description  cf  the 
vorti cl ty-di vergence  shallow  water  model.  Here  the  equations 
are  shown  and  written  using  the  Salerkin  method.  A 
discussion  of  the  computational  technique  used  is  presented 
along  with  the  model's  physical  parameters. 

Chapter  IV  presents  a  descriptive  overview  of  the 
computer  implementation.  The  chapter  includes  a  list  of 
options  available  for  testing,  a  brief  description  of  the 
matrix  compaction  technique  and  the  formulations  of  the 
versatile  modules  used  to  Implement  the  complex  equations. 

Chapters  V  through  VII  discuss  the  results  obtained  from 
the  different  experiments.  Chapter  V  briefly  describes  the 
primitive  model  used  for  all  comparisons  and  the  results 
from  changing  the  element  shape  to  equilateral  triangles. 
Chapter  VI  reformulates  the  primitive  model  so  that  the 
geopotential  is  staggered  with  respect  to  the  velocity 
variable.  7or  simplicity,  the  continuity  equation  is  also 
linearized.  Chapter  VII  compares  the  results  from  the 
vorticity-di vergence  model  developed  in  Chapter  III  to  those 
from  the  primitive  model. 

The  last  chapter  summarizes  the  results  from  all  the 
experiments  and  identifies  what  areas  require  follow  on 
work.  The  source  code  for  the  vorticlty-divergence  model  is 
presented  in  Appendix  A. 
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II.  FINITE  ELEMENTS 

As  is  often  the  case  with  an  original  development,  it  is 
rather  difficult  to  quote  an  exact  date  on  which  the  finite 
element  rrethod  was  invented,  but  the  roots  of  the  method  can 
be  traced  back  to  these  separate  groups:  applied 
mathematicians,  physicists  and  engineers.  Since  the  early 
developments  of  the  finite  element  method,  a  large  amount  of 
research  has  been  devoted  to  the  technique.  However,  the 
finite  element  method  obtained  its  real  impetus  by  the 
independent  developments  carried  out  by  engineers.  Its 
essential  simplicity  in  both  physical  interpretation  and 
mathematical  form  has  undoubtedly  been  as  much  behind  Its 
popularity  as  is  the  digital  computer  which  today  permits  a 
realistic  solution  of  even  the  most  complex  situations. 

The  name  "  finite  element  "  was  coined  in  a  paper  by 
P.V.  Clough,  in  which  the  technique  was  presented  for  plane 
stress  analysis,  as  discussed  in  Bathe  (1976).  While  finite 
element  methods  have  made  a  deep  impact  via  the  field  of 
solid  mechanics,  where  it  can  be  said  that  today  they 
represent  the  generally  accented  method  of  discretizing 
continuum  problems  for  computer-based  solution,  the  same 
appears  not  to  be  true  in  fluid  mechanics  or  atmospheric 
pred iction . 
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Numerous  finite  element  formulations  are  currently 
available.  Strang  (1973),  Norrie  (1373)  and  Zienkiewicz 
'1971)  present  detailed  theoretical  discussions  of  each.  The 
Galerkir  method,  the  most  popular  finite  element  method,  is 
described  in  detail  below  and  used  in  the  equation 
formulation  later. 

A.  FINITE  ELEMENT  CONCEPT 

The  problem  of  solving  partial  differential  equations 
can  be  specified  in  one  of  two  ways.  In  the  first,  finite 
difference  methods  specify  the  dependent  variables  at 
certain  grid  points  in  space  and  time,  and  the  derivatives 
are  evaluated  using  Taylor  series  approximations.  Secondly, 
the  calculus  of  variation  requires  the  minimization  of  a 
functional  over  a  domain,  where  a  functional  is  defined  as  a 
variational  integral  over  the  domain. 

The  calculus  of  variation  approach  creates  a  purely 
physical  model  where  the  functional  equivalent  to.  .the  known 
differential  equations  are  known.  Its  major  disadvantage  is 
that  It  limits  the  method  only  to  those  problems  for  which 
functionals  exist.  Finite  element  methods,  an  extention  of 
this  method,  derive  mathematical  approximations  directly 
fror  the  differential  equations  governing  the  problem.  The 
advantage  here  is  that  it  extends  the  method  to  a  range  of 
problems  for  which  a  functional  may  not  exist,  or  has  not 
been  discovered. 
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The  finite  element  method  divides  the  domain  Into 
subdomains  or  finite  elements  (usually  of  the  same  form). 
Nodes  are  located  along  the  boundary  of  the  elements, 
usually  at  the  element  vertices  and  at  strategic  positions 
(midside,  centroid,  etc.)  in  the  interior  and  on  the  sides 
of  faces  of  elements. 

Commonly  used  elements  are  triangular,  polygonal  or 
polyhedral  in  form  for  two-dimensional  problems.  The  choice 
o'  elements  depends  on  the  type  of  problem,  the  number  of 
elements  desired,  the  accuracy  required  and  the  available 
computing  time.  To  begin  with,  the  element  must  be  able  to 
represent  derivatives  of  up  to  the  order  required  in  the 
solution  procedure,  and  to  guarantee  continuous  first 
derivatives  across  the  element  boundaries  to  avoid 
singularities.  Triangular  elements  are  employed  in  this 
thesis  because  they  can  be  used  effectively  to  represent 
irregular  boundaries,  and/or  geometry,  and  also  to 
concentrate  coordinate  functions  in  those  regions  of  the 
domain  where  rapidly  varying  solutions  are  anticipated. 

Consider  the  problem  of  solving  approximately  the 
differential  equation 

L(u)  =  f (x )  11-1 

where  L  is  a  differential  operator,  u  the  dependent 
variable,  and  f(x)  is  a  specified  forcing  function.  Suppose 
that  I 1—1  is  to  be  solved  in  the  domain  a  *  x  *  b  and  that 
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appropriate  boundary  conditions  are  provided.  The  residual  ?. 
is  forced  from  I I — 1  as  follows: 


L(u)  -  f(x)  *  8 


1-2 


The  critical  step  is  to  select  a  trial  family  cf 
approximate  solutions  (the  members  of  a  trial  family  are 
often  called  basis  functions).  The  basis  function  is 
prescribed  'functionally )  over  the  domain  in  a  piecewise 
fashion,  element  by  element,  and  are  generally  a  combination 
of  low  order  polynomials.  A  one  dimensional  example  is 
shown  in  Figure  1,  wherein  the  domain  (x  axis)  is  divided 
into  six  elements  (line  segments)  A  through  F.  The  basis 
functions  are  linear  and  one  is  shown  for  node  4  only  in 
Figure  1.  The  function  has  a  value  of  unity  over  node  4,  and 
decreases  linearly  to  is  zero  at  nodes  3  and  5  and  zero 
elsewhere. 

Consider  a  series  of  linearly  independent  basis 
functions  Vj  (x),  as  in  Figure  1.  Now  u(x)  can  be 
approximated  with  a  finite  series  as  follows: 


u(x) 


J  ’  Vj 


II  -3 


where  5^  is  the  coefficient  of  the  Jth  basis  function  and  has 
a  value  equal  to  u  at  node  J. 

Substituting  this  approximate  solution  II-3  wherever  u 
appears  in  the  differential  equation  I I — 1 


-  f(x)  -  a 


II -4 
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The  best  solution  will  he  one  which  in  some  sense 
reduces  the  residual  R  to  a  minimum  value  at  all  points  in 
the  domain,  ry  definition,  the  residual  obtained  using  the 
exact  differential  equation  is  identically  zero  everywhere. 
The  residual  P. ,  formed  in  equation  II-4,  is  minimizied  when 
multiplied  with  a  weighting  function,  integrated  over  the 
domain  and  set  equal  to  zero.  This  process  is  known  as  the 
weighted  residual  method 

t 

|  RW  dx  *  3  Il-fi 

a 

where  W  is  the  weighting  function  and  is  referred  to  as  the 
test  function  in  the  following  development.  The  weighted 
residual  method  minimizes  the  errors  of  the  residuals,  such 
that  the  summation  of  all  the  positive  and  negative  errors 
add  to  zero. 

The  Galerkin  method,,  the  most  popular  finite  .gleme at 
method,  is  more  general  in  application  and  is  a  special  case 
of  the  method  of  weighted  residuals,  as  discussed  I/7  Pinder 
and  Gray  (1977).  The  requirement  imposed  on  the  weighted 
residual  method  forming  the  Galerkin  method  is: 

*  the  test  (weighting)  function  he  equal  to  the 
basis  (trial)  function  W  =  V  .  This  process 
leads  in  general  to  the  best  approximation  of 
the  solution. 
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The  final  Galerkin  form  is  obtained  by  substituting  II-4 
into  11*5,  yielding 

b  b 

J  )dx  -  Jvit(z)iz  *  0  1 1 -6 

a  a 

If  this  procedure  is  repeated  for  N  points  in  the  domain 
a  system  with  N  equations  and  N  unknowns  will  be  generated. 

3.  GALEPKIN  application 

The  following  example  taken  in  part  from  Haltiner  and 
Williams  ( 1S80)  applies  the  Galerkin  method  to  the  advection 
equation  with  linear  elements 


du  du 

—  c  — 

at  ax 

=  z 

II-7 

This  equation 

is  dependent 

in 

both  time  and 

space. 

The 

treatment  of  time 

variation 

is 

important 

for 

most 

meteorological  prediction  problems.  The  Galerkin  method  is 
not  applied  to  the  time  dependence  because  it  is  more 
convenient  to  use  finite  differences  in  time, as  is  done  with 
this  example  later.  The  same  treatment  is  applied  tc  the 
prognostic  equations  later,  where  two  finite  differencing 
methods  are  employed  to  do  the  time  integration. 

The  Galerkin  procedure  represents  the  dependent  variable 
ufx.t)  with  a  sum  of  functions  that  have  the  prescribed 
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spatial  structure  as  in  Figure  1.  Approximate  u(x,t)  with 
the  finite  series  as  follows 


u  ( x ,  t ) 


i  4 At )v  (X) 

j=i  J  J 


1 1 -8 


where  the  coefficient  0  .  (t),  a  function  of  time,  is  the 
scalar  value  of  u  at  node  j.  The  basis  functions,  V.(x),  are 

J 

functions  of  space  only  and  J  equals  1  to  7  for  the  example 
in  Figure  1.  The  repeated  subscript  in  this  form  implies  a 
sum  ever  the  repeated  subscript. 

The  Galerkin  equation  for  the  advection  equation  II-7  is 
obtained  by  setting  L  =  c(d()/8x)  and  substituting  in  the 
approximate  solution  II- 8  wherever  u  is  found. 


t  b 

K  (  N  r  8V. 

I  — J  V.V.  dx  +  c  1  0,  — ’ 3  7,  =  0  1 1 -9 

j-i  at  ;  J  1  j-i  J  J  ax  1 

a  a 


where  i  *  1  to  N,  7,  the  test  function  and  V.  the  basis 

^  J 

function.  The  domain  of  integration  is  given  by  a  *  x  *  b, 
and  the  integration  is  done  in  a  piecewise  fashion,  element 
by  element. 

In  this  one-dimensional  case,  an  equation  like  II-S  is 
written  for  each  node,  i.  Considering  node  4,  what  are  the 
possible  non-zero  contributions  from  equation  I 1—97  Figure 
2  illustrates  the  basis  and  the  test  function  interaction 
during  the  piecewise  integration  process.  From  the 
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definition  cf  the  basis  and  the  test  function,  locally 
defined  as  unity  at  node  J  and  linearly  decreasing  to  zero 
at  ?  i  1  and  zero  elsewhere,  the  only  non— zero 
contributions  are  irade  when  *  =  3  over  element  C,  j  *  4  over 
elements  ?  and  r,  and  J  *  5  over  element  3. 

The  evaluation  of  I 1-9  for  i  =  m,  which  is  given  in 
Ealtirer  and  Williams  (1981),  leads  to  the  equation: 


The  boundary  points,  which  in  this  example  are  nodes  1 
and  7,  are  evaluated  in  the  same  way  as  the  interior  nodes, 
with  the  exception  that  cyclic  conditions  are  imposed. 

The  time  discretisation- of  11-10  is  done  using- a—  finite 
difference  scheme.  Applying  leapfrog  time  differencing  gives 
the  following  equation 


The  resultant  equation  set,  in  matrix  form,  contains  an 
NxN  matrix  where  N  is  the  number  of  nodes. 

The  transition  from  one-dimension  to  two  is 
mathematically  identical.  The  domain  is  now  subdivided  into 
finite  areas,  which  are  triangles  in  this  implementation  and 
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the  basis  functions  are  linear.  However,  now  they  are 
pyramid  shaped  with  value  unity  at  the  center  and  decrease 
to  zero  at  the  surrounding  nodes,  and  are  zero  elsewhere. 
Figure  3  shows  this  basis  function  for  node  28  outlined  in 
heavy  black.  The  value  at  any  node  again  can  be  approximated 
by  II-3,  where  J  ranges  over  all  nodes  connected  to  node  i 
including  i  itself.  The  connectivity  for  node  i  *  28  in 
Figure  2  is  J  =  15,16,27,28,29,39  and  43. 

The  integration  is  still  over  the  entire  domain.  With 
both  the  basis  and  the  test  function  zero  over  the  domain, 
except  locally  over  each  element,  the  global  integration  can 
be  performed  by  integrating  locally  over  each  element.  By 
definition,  this  integration  can  be  expressed  as  an  inner 
product  of  both  functions  (i.e.  basis,  test)  as  follows: 

A 

Using  this  definition  and  the  repeated  subscript 
notation  equation  II-9  becomes 

where  the  dot  implies  differentiation  with  respect  to  time, 
and  the  second  subscript  implies  differentiation  with 
respect  to  the  second  subscript.  The  local  integration  may 
be  calculated  directly  from  exact  expressions  derived  from 
area  coordinates  described  in  detail  in  the  next  subsection. 


;/■ 


VidA 


11-12 


29 


Figure  3*  Basis  function  for  node  28.  The 
shaded  area  is  the  coaplete  basis  function 
and  the  Vj  ,  where  j  -  15.  16,27.28,29.39.^ 
are  jth  node  basis  functions  for  node  28.  The 
dashed  line  at  node  28  has  length  unity. 
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In  summary,  the  Galerkin  procedure  Involves  subdividing 
the  domain  Into  finite  elements,  approximating  the  dependent 
variables  by  a  linear  combination  of  low  order  polynomials 
and  substituting  them  into  the  equations.  The  equation  is 
multiplied  by  a  test  function.  Integrated  over  the  entire 
domain  and  finally  the  resulting  system  of  equations  Is 
solved . 

C.  AREA  COOREINATES 

While  the  Cartesian  coordinate  system  Is  the  natural 
choice  of  coordinates  for  most  two  dimensional  problems,  It 
Is  not  convenient  when  working  with  triangularly  shaped 
elements.  It  Is  therefore  necessary  to  define  a  special  set 
of  normalized  coordinates  for  a  triangle.  Area,  or  natural 
coordinates  as  they  are  commonly  called,  reduce  the 
formidable  task  of  Integrating  products  between  the  basis 
and  test  functions  and  their  derivatives  over  a  triangular 
element  and  result  In  easily  computable  -and  exact 
expressions. 

The  following  development  is  taken  in  part  from  the 
formulation  by  Zlenklewlcs  (1971).  Consider  the  triangular 
element  Illustrated  in  Figure  4.  There  is  a  one-to-one 
correspondence  between  the  Cartesian  coordinates  (X,Y)  and 
the  area  coordinates  (L^.Lg.L^  )  for  the  element.  Let  A 
denote  the  area  of  the  triangle  and  A^  ,  A,  and  A^  the  areas  of 
the  subtriangles  In  Figure  4  such  that  A  * 
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The  relationship  between  a  point  P(X,D  in  Cartesian 


coordinates 

and  P(L1# 

l2,  lj)  in 

area 

coordinates 

can  be 

seer,  by  the 

following 

transformations 

X  *  L 

xx  +  l2x  * 

LjX 

Y  *  L 

iY  *  + 

V 

11-14 

1  -  L 

1  *  L2  +  L3 

where 

Li  “■ 

rl  •  l*- 

^  2 

—  and 

A 

l  »i3 

L3  A 

and 

L1  * 

'2k  ♦  btX  ♦ 

ail) 

/2k 

L2  = 

(21  ♦  b2I  + 

a2Y) 

/2k 

11-15 

L3  = 

<2k  *  b3X  - 

a3T) 

/2k 

where  2a  is  twice  the  area  of  the  triangle  and  the  a's  and 
b's  are  defined  as  in  Figure  5. 

It  is  worth  noting  that  every  tuple  ( L  ^  •  L2  ,  ) 

corresponds  to  a  unique  pair  (X,Y)  of  Cartesian  coordinates. 

In  Figure  4,  1^*  1  at  vertex  1  and  0  at  vertices  2  and  3.  A 
linear  relation  exists  between  the  area  and  Cartesian 
coordinates  which  implies  that  values  for  I^vary  linearly 
over  the  triangle  with  a  value  one  at  vertex  1  and  a  value 
of  zero  at  vertices  2  and  3?  and  similarly  for  L2and  L3* 
This  demonstrates  how  each  component  in  the  tuple  (L^,  L2,  L-j) 
behaves  over  the  triangle  as  do  the  linear  basis  and  test 
functions  over  the  element,  as  was  seen  in  Figure  4.  Clearly 

Lt  *  Vi  11-16 
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As  explained  earlier  (see  Equation  11-16),  V^s  a  linear 
function  'l.e.  basis,  test)  which  equals  a  component  Lj.  of 
the  area  coordinate  tuple.  Therefore 


Consequently 


r 


tv. 

1  0 

if  l  A  j 

— J  »  < 

11-21 

bL1 

L  1 

If  1  =  j 

*7j 

for  J  3  1  is 

11-22 


As  an  example,  consider  the  Inner  product  at 

vertices  J  =  2,  1  *  1.  This  Integration  Is  evaluated  as 


A  11-23 

b,  1!  0!  0!  b, 

*  -»  -  2A  *  - 

2A  (1+0  +  0+  2)!  6 

Therefore  any  Inner  product  In  the  formulation  can  be 
readily  evaluated  using  area  coordinates.  Another  benefit  of 
using  this  coordinate  system  is  that  all  of  the  inner 
products  are  functions  of  space  only  and  need  be  computed 
only  once. 
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III.  SHALLOW  WATER  mCISI 

The  governing  equations  for  this  model  are  derived  by 
making  several  simplifying  assumptions  on  the  primitive 
equations  of  motion,  whi oh  then  give  the  barotropic  shallow 
water  equations.  Eowever,  as  mentioned  previously,  the 
shallow  water  equations  describe  many  significant  features 
of  the  large-scale  motion  of  the  atmosphere,  and  therefore 
have  been  used  in  numerous  experiments  over  the  years. 

The  vortlclty-dlvergence  form  of  the  equations  has 
several  advantages.  Williams  (1961)  has  shown  that  the 
geostrophic  adjustment  process  is  treated  much  better  with 
the  vcrticity  divergence  formulation  than  with  a  direct 
treatment  of  the  primitive  form  of  the  shallow  water 
equations,  such  as  was  used  by  Kelley  (1976).  This 
formulation  also  allows  the  velocity  components  and  the 
height  to  be  carried  at  the  same  nodal  points,  whereas  the 
best  scheme  for  the  primitive  form  of  the  equations  requires 
staggering  of  the  fields,  as  seen  in  Schoenstadt  (1980).  The 
vorticity  divergence  form  of  the  equations  is  also 
convenient  for  the  application  of  semi-implicit 
differencing,  which  saves  considerable  computer  time. 
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A.  GOVERNING  ECUATIONS 


The  primitive  form  of  the  shallow  water  equations  in 


Cartesian  coordinates  is 


30  3  ,  *  . 

—  +  73$  *  -  — { 0u)  -  — (0v) 

at  3x  3y 

3u  30  3K 

dt  dx  dx 

3v  30  3K 

dt  djr  &y 


1 1 1 -1 


I II -2 


1 1 1-3 


Equation  f III— 1 )  is  the  continuity  equation  and  the  III-2 
and  II 1-3  are  the  momentum  equations,  respectively.  The 
variables  are  defined  as  follows: 

x ,j  -  the  spatial  coordinates  of  the  domain 

u,v  -  components  of  the  wind  vector 

0  geopotential  »  (gravity  x  free  surface  height) 

2  2 

$  -  mean  geopotential  *  49,000  meters  /seconds 

t  time 

5  -  kinetic  energy 

Q  -  absolute  vorticity  ■  (3  +  fj 

3  relative  vorticity 

f,  -  coriolis  force  (mid-channel  f-plane) 

E  divergence 

The  shallow  water  equations  can  be  written  in 
vorticity  divergence  form  as  follows: 


1 1 1  -4 


a 

a  ; 

—  +  D<p  = 

-  —  (0u) 

- (0v) 

1 1 1  -4 

M 

Bx 

By 

as 

a 

a 

—  s 

- (uQ) 

- (vC) 

1 1 1-5 

at 

Bx 

By 

BE 

a 

a 

r—  *  • 

- (vC) 

- (uQ)  -  v2S 

I II -6 

at 

ax 

By 

where  III-4  is  the  sere  continuity  equation  as  III-l,  III —5 
is  the  vorticity  equation  and  1 1 1 —  €  is  the  divergence 
equation . 

3ecause  of  the  vorticity  divergence  forrr  of  the 
equations,  it  becomes  necessary  to  solve  the  time  dependent 
variables  5  and  I  in  terms  of  ¥,  the  stream  function 
(rotational  part  of  the  wind),  and  1,  the  velocity  potential 
(divergent  part  of  the  wind).  The  initial  fields  for  the 
model  will  be  in  terms  of  %  and 

The  following  diagnostic  relationships  are  defined  and 
used  later  in  the  solution  of  the  equation  set. 


u  *  -  ify+  III-? 
v  =  Sy,  1 1 1-8 
where  the  subscript  implies  differentiation. 


I  *  - 1  "•  hinetic  energy.  III  —9 

2 

uC  *  uOS  ♦  fj,  111-10 

vQ  *  v(*9  +  fj,  III-ll 

«  *  0u,  III  12 

P  -  0v,  I I 1-13 
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*  =  v“  'V, 


1 1 1-14 


D  *  7  *  1, 


I I 1-1= 


3.  2CUAT  TON  FOP.rULATION 

The  Salerkir.  method  described  in  Chapter  II  is  now 
applied  to  eauatlons  III-4  through  III-15.  For  ease  of 
corpr ehension .  the  shorthand  inner  product  dotation  as  in 
11-12  will  he  used  to  simplify  the  equations.  The  detailed 
Galerkin  formulation  will  be  shown  for  equation  III-?,  the  u 
component  of  motion.  The  method  follows  directly  from  the 
example  in  Chapter  II  of  this  thesis,  which  in  turn  follows 
in  part  from  Kelley  (1S76)  and  Ealtiner  and  Williams  (1981). 

Consider  equation  III-?  and  assume  that  each  variable  u, 
^  and  is  auproximated  by 


u  *  V 

*  *  v 

vj  • 


m-ie 


where  the  repeated  subscripts  indicate  summation  over  the 
range  of  the  subscript.  Substituting  these  approximate 
solutions  into  I I 1—7  yields 


d  d 

u,  7  *  —V.  V.)  ♦  — 

jj  dyJJ  9x  J  ^ 


I I 1-17 


Since  only  the  basis  function  is  a  function  of  space, 
III-l?  may  be  further  simplified  by  factoring  out  the  time 
dependent  coefficients. 


39 


The  next  step  requires  multiplying  by  a  test  function  7^ 
as  discussed  in  Chapter  II,  and  integrating  over  the  area 
doma in 

“iff  Vi4A  *  -  *jff  di 

&  A 

*  If  £'i  “ 

A 

The  final  form  in  inner  product  notation  is 

VJ'tl>  -  -  «j  7jy  •5i>  *  <’Cj7j*-Vi>  m*1S 

where  the  double  subscript  implies  differentiating  with 
respect  to  the  second  subscript. 

The  three  prognostic  equations  (III-4,  I I 1—5  and  III —6 ) 
are  similarly  advanced  using  the  Galerkin  technique  to 
become,  respectively: 

*  4  <ejij-,i>  -  "  <Vj*'V  '  ^jW  !I!’szi 

<Vj-V  *  -<(u9)jV’ i>  '  <(rt)jV'i>ni'21 

<VrV  '  ■  <(’s)j  W  - 

+  <Ejv27j,7i>  II 1-22 

2 

where  v  is  the  Laplacian  operator  and  the  dot  implies 
differentiation  with  respect  to  the  time  dependence  in 
III-4,  III- 5  and  III-6. 

Similarly,  Galerkin  equations  are  formulated  for 
equations  III-7  through  III-15. 
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C.  TIKI  EISCRETIZATION 

The  equation  set  III-20,  1 1 1—21  and  III-22  is  arranged 
so  that  all  the  terms  on  the  left  hand  side  nan  be  treated 
implicitly,  and  all  the  terms  on  the  right  hand  side  can  be 
treated  explicitly.  The  explicit  time  integration  will  be 
done  by  the  leapfrog  difference  method.  To  start  the  time 
integration,  two  forward  half  steps  are  taien,  after  which 
the  full  leapfrog  scheme  is  used  for  the  remainder  of  the 
forecast  period. 

The  vorticity  equation  III- 21  is  solved  Independently 
from  I I 1-22  or  III-22.  However,  I 11-20  and  I I 1-22 
'continuity  and  divergence  equations,  respectively)  are 
coupled.  To  explicitly  solve  either,  decoupling  of  the 
equations  is  necessary.  In  this  thesis  this  is  done  through 
algebraic  substitution  of  III-22  (solved  for  D(n+1))  into 
IIT-20.  Once  the  time  integration  is  performed  on  III-20, 
III-22  can  be  solved  for  D(n+1)  using  the  i  (n+l)  value. 

The  final  prediction  equations  are 


t*-' lx>  *  <*.(,. *t,>  *  - 

-  [BEHI]n+1  -  [BERYj  11-1 

*  =^»j  -v  -«rVi> 

-  rfj'k’  j*'vix>  *  <7 jy'vir>’ 

'  -  stlTCjJ^Tp  -  <«c>3<VTi 

-  *C(0u>5<¥3x.Tl  >  *  >} 

I 11-23 
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where  A  =  4/'2*t)  ,  B  *  A/$  ,  C  =  B/(2dt)  and  [BERT]  is  the 
geostrophic  boundary  contribution,  see  Section  3. 

"2j+1<Vj,Vi>  =  s""1  <"j  •^i> 

-  2At  KuC  )  >  *  (vQ)nj<V.yl71  >]  1 1 1-24 

rfS’V  *  Enj~1<?j’vi>  +  (W2)[0n/WJvjfvi> 

-  «°;1<,2'rj'vi>  ♦  2^«>°<vvi> 

-  2  ( uQ  )  j  <Vjy*  V  i  >  +  •2l"<72Vj,V1>J  1 1 1-25 


After  these  three  elliptic  equations  are  solved,  the 
history  of  the  variables  III-?  through  I I 1—15  is  updated. 

A  large  time  step  can  be  applied  to  this  form  of  the 
shallow  water  equations  due  to  the  semi-implicit  nature  of 
the  equations,  -his  is  very  important  since  finite  element 
methods  generally  require  more  computer  time  per  time  step. 
The  vcrticlty-divergence  formulation  acts  as  a  filter,  which 
slows  down  the  high  frequency  waves  in  the  solution.  The 
two-dimensional  advectlve  stability  criterion  for  a  linear 
element,  derived  by  Cullen  (1973),  was  used  to  determine  the 
correct  time  step. 


4  X 

Icl  fe 


I I 1-26 


whereat  is  the  time  step  in  seconds,  *x  the  shortest  grid 
spacing  in  meters  and  c  the  fastest  phase  velocity. 


E .  COMPUTATIONAL  TECHNIQUES 

The  final  prognostic  equation  set  requires  the  solution 
of  a  Helmholtz  equation  for  $  and  Poisson  equations  for  ¥ 
and  7.  The  irost  common  method  of  solution  used  by 
rreteorologists  has  been  the  successive  over  relaxation 
method  (SOR)  in  which  an  Initial  guess  of  the  solution  is 
made  and  then  progressively  improved  until  an  acceptable 
level  of  accuracy  is  reached.  SOR  is  employed  in  the 
solution  of  the  equations, where  I I 1-23  can  be  represented  by 

v2  [M]  { x}  -  C[M]{x}  -  {b}  III-27 

and  III  24,  III  25  by 

v2[M]{x}  *  {b}  II 1-28 

where  v2  the  Laplacian  operator,  [M]  *  <V  j»7i>  Matrix,  {x} 
-  the  dependent  variable  in  vector  notation,  C  -  constant  as 
in  III -23  and  {b}  the  right  hand  side  of  the  equation  or  the 
forcing  function. 

The  mass  matrix  [M] ,  dimensioned  (nxn),  is  a  matrix  of 
coefficients  whose  rows  are  the  equations  of  the  system  to 
be  solved.  There  exists  a  one  to  one  correspondence  between 
the  rows  of  the  mass  matrix  and  the  nodes  of  the  domain. 
Each  equation  has  a  term  (column)  for  each  node,  where  a 
non-zero  term  represents  connectivity.  Nodes  are  connected 
if  they  are  both  vertices  of  the  same  element.  Obviously  CM] 
is  a  sparse  matrix  containing  the  inner  products  for  the 
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left  hand  side.  Chapter  4  of  this  thesis  will  describe  the 
matrix  compaction  procedure. 

The  forcing  function  {t},  dimensioned  (nxl),  involves 
only  variables  at  the  current  time  step  and  is  easily 
computed  using  four  very  versatile  subroutines  described  in 
detail  in  the  next  chapter. 

The  initial  guess  to  start  SOS  is  the  previous  time  step 
solution.  An  average  of  30  passes  per  equation  are  needed 
for  each  time  step.  The  solution  is  considered  to  have 
converged  to  its  final  value  when  the  residual  for  each  node 
has  been  reduced  to  some  acceptably  small  value. 

The  diagnostic  equations  III-7  through  III-15  must  also 
be  solved  every  time  step.  However,  the  same  technique  is 
not  used  for  these  equations.  Dr.  **.J.P.  Cullen  suggested  an 
under  relaxation  scheme  for  which  three  passes  over  the 
domain  should  produce  a  solution  of  acceptable  accuracy, 
since  the  coefficient  matrix  is  so  strongly  diagonally 
dominant,  Mass  lumping  cf  the  coefficient  matrix  is  used  for 
the  first  guess.  This  technique  requires  replacing  the  mass 
matrix  M  by  the  identity  matrix  [I].  A  first  guess  of  this 
type  is  able  to  describe  most  of  the  large  scale  features, 
which  in  turn  reduces  the  number  of  iterative  passes  over 
the  field.  Successive  passes  converge  to  solutions  which 
describe  smaller  scale  motion,  approximately  to  the  same 
order  of  magnitude  as  introduced  by  computational  errcr,  so 
that  further  iterations  are  not  needed. 
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2.  GRID  GEOMITRT 

The  domain  of  this  model  is  a  cylindrical  channel,  with 
total  length  of  4045  Km  and  width  of  3503  Km.  The  channel 
simulates  a  belt  around  the  earth  and  it  proves  to  be  an 
excellent  test  bed  for  comparing  with  the  finite  element 
formulations  used  by  Kelley  (1576)  and  Older  (1981). 

The  domain  is  subdivided  into  equilateral  triangles  as 
shown  in  Eigure  6.  most  of  the  test  runs  for  this  thesis  use 
a  12x12  mesh  which  has  156  nodes  and  289  elements.  This 
implementation  is  not  restricted  to  one  grid  pattern.  The 
technique  developed  by  Older  (1981)  to  vary  the  nodal 
geometry  smoothly  to  achieve  areas  of  denser  and  coarser 
resolution  is  also  implemented,  as  in  a  third  grid  pattern 
that  varies  the  nodal  geometry  abruptly.  A  short  discussion 
of  these  nodal  geometries  with  accompaning  illustrations  of 
each  is  presented  in  Chapter  VII,  where  the  different  test 
cases  are  described. 

Cyclic  continuity  is  assured  in  the  x  direction  by 
wrapping  the  domain  around  the  earth  to  form  a  cylindrical 
domain.  This  has  the  advantage  of  eliminating  the  east-west 
boundaries  and  it  simulates  the  flow  around  the  earth.  The 
only  boundaries  c n  this  domain  are  the  north-south  walls  and 
their  treatment  will  be  discussed  shortly. 
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Figure  6  Domain  divided  into  equilateral  triangles. 


F.  INITIAL  CONDITIONS 


As  mentioned  previously,  the  reformulation  of  the 
frcverning  equations  into  the  vort  icity-di7ergence 
shallow  water  equation  set  requires  solving  the  time 
dependent  variables  in  terms  of  the  stream  function  and 
velocity  potential.  The  continuity  equation  is  not  altered, 
so  that  its  solution  is  expressed  in  terms  of  0. 

For  the  basic  testing  of  the  model's  performance,  simple 
analytic  sinusoidal  initial  conditions  are  used  to  insure 
the  rost  accurate  analysis  possible  and  to  simplify  the 
comparisons. 

The  sinusoidal  initial  fields  are  graphically  shown  in 
Figure  7  as  3-dimensional  surfaces.  The  geopotential  field  0 
consist  of  a  half  sine  wave  in  the  y  direction  and  a  single 
cosine  wave  in  the  x  direction.  The  stream  function 
calculated  by  dividing  the  geopotential  field  by  the 
coriolis  force,  has  the  same  physical  structure  as  0.  The 
velocity  potential  X  has  a  .  single  sine  wave  in  the  x 
direction  and  a  half  sine  wave  in  the  y  direction. 

These  initial  conditions  are  computed  as  follows 


0 

=  fo  Asinoc^cosoCg  ~ 

■  f„a<jr  -  Vj  -  i 

¥ 

=  0/fo 

III -29 

X 

*  CsinocjSinafg 

quasi-geostrophic 

divergence 

where 

A  -  arbitrary 

amplitude 

fc  -  coriolis 

value  for  mid-channel  latitude 

a)  and  ¥  initial  fields. 


b)  X  initial  field. 


Figure  .  3_dio«nsional  view  of  the  inital  fields. 


U  mean  flow 
ya  -  mid-latitude  value  of  y 
|  -  mean  free  geopotential  height 

=  49,200  m2/s  2 

otj_  -  iry/V 
o<2  2mxr/L 
r  -  wave  number 
V  channel  width 
L  -  channel  length 
C  -  -(f0Uoc2BA)/(fS  ♦  $2) 

B  -  Or*  +  <*2 

S.  BOUNEAHY  CONEITIONS 

Boundary  conditions  are  only  required  on  the  north  and 
south  walls  of  the  grid  dorrain.  Cue  to  cyclic  continuity, 
the  domain  is  wrapped  around  creating  a  cylinder  eliminating 
the  east  and  west  boundaries.  However,  careful  attention  to 
detail  is  needed  during  the  implementation  to  assure  this 
continuity.  Separate  boundary  conditions  are  applied  to  each 
of  the  predictor  equations  III-23,  III-24  and  III -25 .  These 
conditions  are  computed  for  the  wall  nodes  only  and  are 
applied  during  each  pass  through  the  relaxation  scheme. 

The  vorticity  equation  III-24,  the  most  sensitive  of  the 
predictor  equations  to  solve,  requires  if  on  the  north-south 
boundaries  to  remain  constant  for  the  entire  forecast 
period.  Since  this  equation  is  solved  in  terms  of  if,  the 
Initial  north-south  if  values  are  saved  and  assigned  to  the 
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boundary  points  after  each  pass  through  the  relaxation 
subroutine. 

'"he  proper  boundary  condition  for  the  divergence 
equation  III —25  would  be  dX/dn  =  Z.  However,  for  the  purpose 
of  this  study,  there  is  more  interest  in  the  sinusoidal 
variation  in  the  y  direction  and  not  in  the  region  of  the 
walls.  Therefore  X  =  2  is  appropriate. 

The  continuity  equation  III-23,  the  most  complex 
predictor  equation,  requires  that  there  be  no  mass  flux 
through  the  north-south  walls.  The  geostrophic  boundary 
condition 

—  »  -uf  1 1 1-30 

ay 

is  auplied  to  the  north  south  boundary  nodes  for  the  terms 
[BERT]  in  equation  III-23.  Integrating  the  inner  product 
7 . ,7, ^  by  parts  produces  the  boundary  terms 

J  J  A 

Jf  v2^jVj)Vi  dxdy  *  f/ *  •  5 Vj_  dxdy 

y  x  yx 

*  f(  Cv*(71T(5(;37j)}  -  dxdy 

yx 

*  i  dr  “  M  77i*7^j7j^  ixd7 

yx 

*  [BDRTl  -  ^[<7  jx,Vlx>  *  <7jy.7iy>]  HI-31 

where  n  is  a  unit  vector  normal  to  the  domain  and  dr  is  the 
differential  distance  along  the  path  of  integration  on  the 
perimeter  of  the  domain. 
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The  geostrophic  boundary  condition  I I 1—32  is  substituted 
into  the  contour  integral  in  equation  III —3 1  and  put  into 
lalerkin  forrr,  in  the  same  way  as  in  the  one-d irrenslonal 
advective  equation  in  Chapter  II.  The  resulting  term  is 
derived  as  follows 


-  ySi^VV’i  m'S2 

zauation  III-32  appears  twice  in  the  continuity  equation 
III  23,  for  time  levels  (n+1)  and  (n-1).  All  values  of  u  are 
known  for  time  (n-1),  since  they  are  saved  from  the  previous 
calculations.  However,  u(n+l)  has  not  been  computed.  To 
solve  for  u(n+l),  both¥(n+l)  and  T(n+1)  are  needed,  ^(n+l) 
is  solved  first  from  the  vorticity  equation.  l(n+l)  needs 
0("+l)  as  part  of  its  solution  and  0{n*l)  needs  u(n+l)  in 
its  solution. .  To  avoid  this  problem,  it. is  assumed  that 
T(n-t-l)  has  a  negligible  contribution  to  the  solution  of 
u(n+l)  and  only  * ( n+1 )  is  used. 


IV-  CCrPUTIP.  IMPLANTATION 


The  formulation  and  general  theory  of  the  finite  element 
method  was  presented  in  the  previous  chapters.  The  objective 
ir  this  chapter  is  to  discuss  sore  important  computati onal 
aspects  pertaining  to  the  implementation  of  the  finite 
element  prediction  system. 

The  main  advantage  that  the  finite  element  method  has 
over  other  prediction  techniques  is  its  generality. 
Conceptually,  it  seems  possible  by  using  many  elements,  to 
approximate  virtually  any  surface  with  complex  boundaries 
and  initial  conditions  to  such  a  degree  that  an  accurate 
solution  can  be  obtained.  In  practice,  however,  obvious 
engineering  limitations  arise,  a  most  important  one  being 
the  cost  of  the  computation.  As  the  number  of  elements  is 
increased,  a  larger  amount  of  computer  time  is  required  for 
a  forecast.  Turthermore,  the  limitations  of  the  program  and 
the  computer  may  prevent  the  use  of  a  large  number  of 
elements.  These  limitations  may  be  due  to  the  computer  speed 
and  storage  availability,  or  round-off  errors  propagated  in 
the  computations  because  of  finite  precision  arithmetic. 
Also,  the  malfunction  of  a  hardware  component,  if  the 
prediction  is  carried  out  using  many  computer  hours  to 
execute,  can  be  a  serious  problem.  It  is  therefore  desirable 
to  use  efficient  finite  element  programs. 
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The  effectiveness  cf  a  program  depends  essentially  or. 
the  following  factors.  Firstly,  the  use  cf  efficient  finite 
element  techniques  is  important.  Secondly,  efficient 
programming  methods  and  sophisticated  use  cf  the  available 
computer  hardware  and  software  are  important.  The  third  very 
important  aspect  in  the  development  cf  a  finite  element 
nrogram  is  the  use  of  appropriate  .numerical  techniques. 

The  vortici ty-divergence  model  described  in  the  previous 
chapter  is  implemented  on  the  131*  3033  computer  located  at 
the  Naval  Postgraduate  School.  Some  notatle  features  of  its 
architecture  are  the  three  trillion  bytes  of  virtual  mass 
storage,  of  which  eight  mega  bytes  are  available  to  each 
user,  and  the  £7  nanosecond  machine  cycle  time.  The  model  is 
executed  mostly  using  a  12x12  element  domain  requiring  4 00k 
bytes  of  storage  and  30  seconds  of  CPU  time  to  execute. 
Exceeding  execution  time  and/or  available  storage  is  not  a 
problem,  in  fact  the  system  allowed  a  lot  of  flexibility 
during  the  implementation  phase  of  the  model. 

The  source  code  is  written  using  FORTRAN  IV  and  compiled 
on  an  optimizing  Fortran  E  compiler.  Appendix  A  contains  the 
source  cede  listing,  which  is  divided  into  five  subdivisions 
delineating  the  logical  structure  of  the  program. 


A.  PROGRAM  ARCHITECTURE 

Program  features  incorporated  in  the  model  are: 

1)  vodularity.  With  only  a  few  exceptions,  each 
module  is  limited  to  one  page  of  FORTRAN  code.  This  makes  it 
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easier  to  comprehend  the  program.  lach  module  performs  only 
one  task.  ?cr  exarple,  subroutine  CONTIC  computes  the  value 
c*’  the  forcing  terms  for  the  continuity  equation.  likewise 
there  is  also  one  module  for  the  divergence  and  vcrticity 
equations.  To  implement  a  new  set  of  equations,  only  these 
red',  les  would  have  to  oe  altered. 

2:  "aslly  controllable  switches.  Switches  may  be  set 
to  either  print,  plot  or  tabulate  harmonic  analysis  data  for 
most  of  the  available  fields.  The  ability  to  display 
intermediate  results  allows  each  portion  of  the  algorithm  to 
be  monitored  for  computational  adjustments.  This  also  makes 
it  easier  for  unfamiliar  users  to  become  acauinated  with  the 
computational  model. 

3  ^  forcing  term  subroutines.  In  previous 
implementations .  each  forcing  term  was  calculated  by  a 
special  subroutine.  In  this  implementation,  the  calculations 
are  accomplished  by  general  purpose  routines  which  simplify 
the  implementation  of  the  complex  prognostic  equations.  Thl-s- 
allows  implementation  of  different  equation  sets  (i.e. 
^arocllnlc  yodel)  over  the  same  domain  with  minimal  effort. 

4)  Documentation.  Each  variable  is  defined  by  a 
short  phrase  (Appendix  A,  A.).  The  function  of  each  module 
is  described  in  an  introductory  paragraph.  Shared  data  is 
placed  in  named  common  blocks  and  identified  with  each 
subroutine  which  uses  them.  A  subroutine  index  is  given. 
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The  main  program  is  short,  calling  only  six  modules 
which  reflect  the  "basic  sequential  flow  cf  the  model.  It 
starts  with  initialization  of  all  model  parameters  (i.e. 
model  options,  domain,  finite  element  arrays,  inner 
products).  It  then  Initializes  the  input  fields  (i.e. 
gecpctential  heights,  stream  function  and  velocity 
potential)  and  is  followed  by  initialization  of  all 
remaining  dependent  variables.  At  this  time  the  model  is 
totally  initialized  and  time  integration  begins.  As 
mentioned  previously  two  techniques  are  employed  for  time 
integration,  each  having  its  own  module.  Upon  completion  of 
the  last  forecast,  the  program  terminates. 

Arrays  are  the  only  data  structures  used  and  are 
grouped  using  IS  different  common  blocks.  Several  arrays  are 
used  as  static  link  lists,  as  described  in  detail  later, 
which  simplified  the  algorithms.  The  common  block  format  has 
the  advantage  of  reducing  the  overall  execution  time  of  the 
program.  Post  of  the  arguments  passed  during  a  call  to  a 
subroutine  are  contained  in  conmon.  This  requires  less  time 
to  execute  since  no  parameter  passing  is  required  for  the 
arguments.  Another  benefit  of  this  format  is  that  the  code 
becomes  less  cumbersome  and  more  readable.  Each  variable  and 
array  Is  defined  in  the  first  subsection  of  Appendix  A  along 
with  a  page  index  for  the  subroutines. 


2.  Initialization  Phase 


Appendix  A,  Section  C  contains  all  the  subroutines 
used  during  the  initialization  phase  of  the  program.  From 
the  user  point  of  view,  the  most  important  suoroutine  is 
the  first  subroutine  called,  which  contains  all  the 
global  variables  that  control  the  different  options 
available  per  run.  This  is  the  only  subroutine  that  is 
changed  to  run  the  different  experiments,  assuming  that  no 
new  computational  technique  is  introduced.  The  selection  of 
options  are: 

1)  -  channel  location  -  the  channel  may  be 

shifted  north  or  south  by  presetting  the 
north/south  latitude  limits  in  INITG3. 

2)  -  variable  geometry  -  the  node  positions  may 

be  grouped  for  more  dense  node  patterns  to 
yield  higher  resolution.  Two  variables  HI 
and  H2  set  the  ratio  used  to  vary  the 
spacing  along  the  x  and  :  y  axes, 
respectively. 

3)  -  initial  field  wave  length  and  amplitude  can 

be  altered  to  produce  various  effects. 

4)  change  the  initial  mean  flow. 

5)  -  diffusion  can  be  entered  for  any  of  the 

three  prognostic  equations. 

6)  -  maximum  length  of  forecast  period  may  be 

changed  and  a  print,  plot  or  harmonic 


55 


analysis  of  any  dependent  variable  may  be 
requested  for  any  time  interval. 

Or.  ce  the  experimen  t  is  dete  rained ,  the  options 
listed  above  are  set.  The  program  is  ready  to  he  executed. 

The  largest  part  of  the  initialization  pnase 
consists  of  establishing  the  domain  and  producing  all  the 
finite  element  computational  vectors  that  remain  constant 
throughout  the  experiment. 

The  first  several  steps  in  setting  up  the  domain  are 
concerned  with  indexing.  Subroutine  COP.RES  is  called  first, 
where  all  the  nodes  (grid  points)  and  elements  (triangular 
areas)  are  numbered  consecutively  starting  at  the  southwest 
corner  of  the  domain  and  moving  eastward  across  each  row  or 
latitude.  For  each  'element;--  a  record  of-ali.  of  its  nodes 
(vertices)  are  stored  in  array  ELK ENT  (T.,3),  where  M  is  the 
total  number  of  elements.  To  facilitate  the  inner  product 
evaluation  later,  a  local  numbering  system  is  required  for 
each  element.  That  is,  for  each  element,  its  nodes  are 
stored  counterclockwise  in  a  positive  sense.  The  first  node 
however,  is  arbitrary. 

Vlth  the  domain  divided  ar.d  numbered,  a  connectivity 
list  (the  correspondence  between  each  node  and  the  neighbor 
nodes)  is  constructed  for  each  node  by  subroutine  CORRIL. 
Each  node  is  adjacent  to  four  or  six  other  nodes  depending 
on  whether  it  is  a  boundary  or  interior  node,  respectively. 
These  adjacent  nodes,  plus  itself,  make  up  the  connectivity 
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list  for  one  node.  The  connectivity  lists  are  then 
conca tena ted  sequentially  starting  with  the  first  nodes 
connectivity  list  into  the  vector  NAME  (NN),  where  NN  is  the 
sur  of  the  nodes  in  each  connectivity  list.  (i.e.  for  a 
12x12  domain  with  156  nodes,  and  equilateral  elements  NN  = 
15*44).  5or  the  first  tire  during  the  initialization  phase, 
special  attention  is  given  to  cyclic  continuity.  As 
discussed  earlier,  cyclic  continuity  is  the  joining  of  the 
east  and  west  boundaries  to  create  a  cylindrical  channel. 
The  connectivity  list  for  these  east/west  boundary  nodes 
rust  be  complete  to  insure  proper  continuity  for  the 
ca  tions  later. 

The  connectivity  vector  NAPE  is  frequently  used 
during  most  computations.  Two  utility  vectors  ISTART 
(containing  the  starting  location  In  NAPE  for  a  particular 
node)  and  MOM  (containing  the  number  of  nodes  in  its 
connectivity  list)  are  used  to  locate  and  index  through  the 
vector  NAME,  as  will  be  seen  shortly.  This  same  technique  is 
used  to  index  through  the  coefficient  matrices  and  used 
during  most  of  the  node  interaction  computations. 

The  physical  properties  of  the  channel  are 
calculated  next  in  subroutine  C5ANAL.  Eere  the  north  and 
south  latitude  boundaries,  which  were  pre-set  in  INITGR  by 
the  user,  are  used  to  compute  the  grid  spacing  along  the  x 
and  y  axis.  Since  this  channel  simulates  a  belt  around  the 
earth,  the  magnitudes  of  both  EELTAX  and  EELTAY  (meters)  are 
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proportional  to  the  width  of  the  channel  divided  hy  the 
number  of  grid  points  in  the  y  aiis. 

The  Cartesian  coordinates  for  each  node  are  computed 
by  subroutine  LOCATE  using  the  DELTA!  and  DELTA!  calculated 
in  CEANAL.  If  the  option  to  use  varying  grid  geometry  is 
desired,  subroutine  TRANS  transforms  the  grid  geometry. 
TRANS  also  computes  the  corresponding  new  Cartesian 
coordinate  values  for  each  grid  point  and  calculates  the 
minimum  IELTaX  and  TELTAY  within  the  domain.  When  the 
geometry  is  changed  to  create  a  smaller  DELTA!  or  DELTA! , 
the  two  dimensional  advective  stability  criteria  is  also 
charged.  A  new  time  step  DT  has  to  be  computed  using 
equation  III-2e.  Since  TRANS  transforms  the  geometry,  it 
also  computes  the  new  DT. 

Another  transformation  is  required  as  discussed  in 
Chapter  II.  The  transformation  from  Cartesian  coordinates  to 
area  coordinates  is  needed  to  perform  the  area  integration 
of  the  inner  products.  Subroutine  AREA  computes  -thes-e 
transformations  exactly  as  outlined  in  Chapter  II,  Section 
C.  Again  cyclic  continuity  is  very  important  and  special 
care  is  needed  to  insure  proper  transformation. 

Following  the  area  transformation  is  the  computation 
of  all  the  inner  products  that  are  required  to  solve  the 
equations.  The  advantage  of  using  area  coordinates  is  that 
the  inner  products  (function  of  space  coordinates  only)  are 
computed  and  stored  once  and  used  repeatedly  without 
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recalculation.  Subroutine  INNER  computes  and  stores  these 
products  using  the  formulas  derived  in  Chapter  II. 

Che  coefficient  matrix,  dimensioned  NxN ,  where  N  is 
the  total  number  of  nodes,  is  a  matrix  of  coefficients  whose 
rows  are  the  equations  of  the  system  to  be  solved.  As 
discussed  in  the  computational  technique  section,  the 
members  of  this  sparse  matrix  are  the  inner  products  for  the 
left  hand  sides  of  the  equations.  Three  coefficient  matrices 
are  used  in  the  solution  of  the  equations.  The  diagnostic 
equations  (III-7  through  III-15)  use  a  coefficient  matrix 
with  the  inner  product  <V  j  ,V^  >  which  is  constructed  by- 
subroutine  AMTRTl  and  stored  in  compacted  form  in  vector 
GfN'N)  by  subroutine  ASE^BL.  However,  when  solving  the 
pro*nostic  equations,  these  coefficient  matrices  have  a  DT 
f time  step  in  seconds)  term,  so  that  these  matrices  are  not 
assembled  until  the  time  integration  begins.  The  vortlcity 
and  divergence  equations  (111-24,111-25)  use  the  coefficient 
matrix  S(NN)  with  inner  products  <7^,7  i 10 
solving  the  Poisson  equations  for  the  stream  function  and 
velocity  potential,  respectively.  The  continuity  equation 
(III-23)  uses  a  combination  of  inner  products  in  its 
coefficient  matrix  F(NN)  as  follows 


<?jx’Vi*>  +  <*jyTiy>  *  ^2<Vj’Vi> 


ea 
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to  solve  the  Helmholtz  equation.  At  the  start  of  each  tire 
integration  module,  subroutine  AMTRI2  is  called  tc  construct 
the  two  coefficient  ratrices  H  and  F. 

These  banded  and  sparse  ratrices  are  ccrpacted  into 
vectors  to  save  storage  during  their  asserblage  by  ASEM3L. 
The  vectors  are  dirensioned  NN,  as  is  NAME,  the  connectivity 
vector,  and  both  use  ISTART  and  MUM  to  index  through  them. 
This  corpaction  routine  was  used  by  Kelley  (1976)  and  Older 
'1931)  in  their  rodels,  but  was  developed  by  Hinsran  (1975). 

To  illustrate  ratrix  asserblage  using  an  elerent  by 
elerent  technique,  consider  Figure  8.  Note  that  this 
illustration  is  for  elerent  number  3,  but  all  elements  are 
treated  in  a  similar  raner.  The  computational  technique 
required  that  for  each  point  (node)  describing  element  3, 
namely  nodes  2,  3  and  14  stored  in  array  ELMENT ,  the  inner 
product  between  those  points  be  distributed  to  their 
proper  location  in  the  coefficient  matrix. 

Subroutine  AMTRX1  builds  the  inner  product  nodal 
Interaction  and  stores  it  in  matrix  3,  dimensioned  3x3. 
Figxre  8  illustrates  the  3  matrix  for  element  3,  where  the 
inner  product  <V  j  »vi>  ls  tlie  multiplicand  of  the 
corresponding  basis  and  test  functions,  respectively. 

The  local  dispensing  of  interactions  is  done  in 
ASEMBL.  Consider  the  second  row  of  [3]  in  Figure  8.  These 
are  the  interactions  between  node  3  of  element  3  to  the  test 
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Figure 


AMTRX  -  builds  [b]  for  one  element,  passes  [2]  and  the  element  number 
associated  with  [3]  to  aSEKBL.  This  process  continues  till  all  element 
node  Interactions  are  assembled  in  the  coefficient  matrix. 
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ASEMBL  -  assemblea  the  node  interactions  into  the  coefficient  matrix 
[c]  ,  which  ham  the  same  structure  as  HAKE.  The  following  diagram  assembles 
inner  product  Into  the  coefficient  matrix  [c] ,  for  element  J, 
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8.  Assembling  and  storing  the  coefficient  matrix  for  element  3. 
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functions.  AS1PBL  locates  nodes  3's  connectivity  list  in 
NAPI  using  ISTART  and  NUP.  In  Figure  S,  this  list  is 
delineated  by  STA?.T(3)  and  LAST(3).  Now  ASEpBL  steps  through 
the  connectivity  list  for  three  iterations.  Turing  each 
pass.  ASZM3L  is  searching  for  one  of  the  three  node  numbers 
for  element  3.  ’#hen  a  iratch  is  found  with  one  of  element  3's 
nodes  'i.e.  2,3  or  14)  and  node  3's  connectivity  list  (i.e. 
3,2,14,15  or  4)  the  proper  position,  to  which  this 
interaction  is  to  be  added  has  been  found  in  the  coefficient 
matrix.  Since  NAPE  and  vector  C,  the  compacted  coefficient 
matrix,  are  dimensioned  identically,  the  same  pointer  (i.e. 
J  in  Figure  e)  is  used  to  index  through  both  arrays.  This 
procedure  is  repeated  for  all  elements  in  the  domain  to 
assemble  the  coefficient  matrix  of  the  equations.  The 
pseudo-code  for  ASEMBL  is  shown  in  Figure  8  to  facilitate 
stepping  through  this  example. 

The  domain  and  all  finite  element  worfc  vectors  are 
initialized  at  .thiSi  patat.  -SuhrautLae  EiP.SJtT.  is  -calted  -Later 
to  compute  interpolation  points  for  the  harmonic  analysis 
subroutines . 

The  last  phase  of  the  initialization  process  is  the 
initialization  of  the  dependent  variables.  The  three  input 
fields  geopotential  heights,  stream  function  and  velocity 
potential  are  computed  in  subroutine  IC  using  the  equation 
set  III-30.  However,  the  variables  calculated  from  the 
diagnostic  equations  have  to  be  computed  using  the  input 
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fields.  These  variables  are  used  during  each  time  step  while 
solving  the  prognostic  equations. 

The  diagnostic  equations  are  solved  in  subroutine 
EEPVAP.,  first  during  the  initialization  phase  and  later 
during  the  time  integration  phase.  Each  diagnostic  equation 
calls  its  own  module  to  compute  the  value  of  the  forcing 
function  and  stores  the  computed  values  in  the  vector  RES. 
These  equations  all  use  the  same  coefficient  matrix  when 
solving  the  diagnostic  equations.  Subroutine  SOLVER  is 
sufficiently  genereal  to  solve  each  equation.  SOLVER  uses 
vector  RES  and  coefficient  matrix  G  to  under-relax  towards 
the  solution.  As  mentioned  previously,  the  coefficient 
matrix  is  strongly  diagonally  dominant  so  that  three  passes 
over  the  domain  are  sufficient.  At  the  end  of  EEFVAR,  output 
is  generated  depending  on  what  print,  plot,  or  harmonic 
analysis  controls  were  preset. 

This  completes  the  initialization  phase  of  the  model 
and  the  program  for  the  forecast  phase  will  be  described 
next . 

3 .  Forecast  Phase 

The  forecast  phase  is  accomplished  in  two  steps.  The 
first  time  step  is  made  using  two  half  steps  by  subroutine 
FATZNO.  Eere  the  prognostic  coefficient  matrices  are 
constructed  using  half  the  BT  value  by  calling  AMTRX2. 
AMTRX2  uses  the  same  computational  technique  to  construct 
the  coefficient  matrices  as  described  for  AtfTRXl . 
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Each  of  the  prognostic  eouations  III-23,  III-24  and 
I 11-25  calls  its  own  subroutine  (CONTEC,  VORTEQ  and  LIVEC 
respectively)  to  compute  all  the  terms  on  the  right  hand 
side,  which  are  stored  in  the  vector  RBS.  After  computations 
for  ?.HS  are  completed,  subroutine  RELAX  solves  the  equations 
by  over-relaxation  as  described  in  the  computational 
technique  in  Chapter  3.  Once  the  solutions  for  the  'n+l) 
time  step  are  completed,  DEPAR  is  called  to  update  the 
variables  from  the  diagnostic  equations.  Two  passes  through 
^ATZNO  advances  the  solution  fields  one  time  step. 

The  remainder  of  the  forecast  period  is  integrated 
using  the  leapfrog  scheme.  Subroutine  LZAPER  performs  this 
Integration  using  the  identical  format  as  MATZNO,  except 
that  E?  equals  two  DT.  At  preset  times  as  specified  in 
INITG2,  the  different  fields  are  saved  for  printing.  This 
process  continues  until  the  final  forecast  time  is  reached. 

s.  utility  moruizs 

Once  the  equation  formulation  is  completed,  as  in 
Chapter  III,  all  the  inner  products  and  types  of 
integrations  are  known.  Versatile  modules  can  be  written  to 
perform  these  computations.  Consider  a  term  of  general  form 
^Aj  7^  ,V  ^  >  where  i  is  the  node  about  which  the  term  is 
evaluated  and  the  j's  are  the  nodes  connected  to  node  i,  or 
the  surrounding  nodes.  The  inner  product  values  are 
already  computed  and  stored  for  all  the  nodes,  during  the 
initialization  phase  of  the  model. 
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"he  regaining  computation  to  complete  tne  evaluation  of 
this  term  is  the  multiplication  of  the  scalar  coefficient  of 
A  at  node  *  with  the  corresponding  inner  product  < 7 . , V, >  for 

s) 

node  ’.  This  requires  indexing  through  node  i 's  connectivity 
list  stored  in  vector  NAME,  and  for  each  node  in  the  list 
multiply  and  total  the  products.  The  cumulative  sum  of  these 
multiplications  is  assigned  as  node  i 's  contribution  for 
this  term.  Subroutine  TERM3  performs  this  exact  computation. 
All  that  is  passed  to  TERM3  is  the  scalar  field  A  and  the 
sign  of  the  inner  product,  TERMS  then  computes  the 
contribution  for  each  node  in  the  domain  and  accumulates  it 
in  the  work  vector  P.HS  . 

Three  other  utility  modules  are»  TIRM1 ,  which  computes 
the  first  scalar  multiplication  for  triple  inner  products 
(i.e.  O.j  Vj  3fcVk  ,V1>) .  The  product  Is  again  already 

computed  and  stored  by  subroutine  INNER.  TE3M1  computes 
vk  .V  i  ^  and  construct*  a  compacted  vector  similar  to  the 
coefficient  matrices.  This  reduces  the  effort  of  multiplying 
the  second  scalar  to  a  TERM3  computation.  TERM2  computes 
node  Interaction  of  the  following  type  <Aj  where  both 

the  basis  and  the  test  functions  are  derivatives.  Lastly, 
subroutine  TIRM4  computes  node  interaction  for  terms  as  <A^ V^x 
, V i .■  f  where  only  the  basis  function  is  a  derivative. 

When  examining  the  right  hand  side  of  the  equation  sets 
III-23,  III-24  and  III-25,  it  is  obvious  this  implementation 
is  a  subscripting  nightmare;  however,  the  use  of  the  utility 
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modules  TESM1 .  TERM2 ,  TERM3  and  TERM4  simplifies  the 
implementation  to  only  determining  what  order  to  call  the 
utility  modules.  Examination  of  subroutines  CONTEQ,  VCRTZQ 
and  EIVEC,  which  compute  the  right  hand  sides  for  III-23, 
111-24  and  1 1 1—25  respectively,  illustrates  this  fact.  No 
other  subroutines  or  calculations  were  required. 
Implementation  of  these  equations  required  minimal  effort. 


7.  PRIMITIVE  MOEZL  EXPERIMENT 


The  previous  two  chapters  presented  the  detailed 
formulation  and  implementation  of  the  vorticity-divergenoe , 
shallow  water  equation  model.  The  results  from  this  model 
will  he  compared  with  the  results  from  the  primitive  model 
in  Chapter  VII.  To  facilitate  interpretation  of  the 
comparisons,  a  brief  description  of  the  primitive  model 
follows.  See  Kelley  (1976)  for  a  detailed  discussion  of  the 
entire  model. 

Also  presented  in  this  chapter  is  an  experiment  which 
demonstrates  significant  Improvement  of  the  solution  from 
the  primitive  model.  Kelley's  implementation  used  elements 
which  were  right  triangles.  Older  (1981)  showed  that 
eauilateral  elements  are  far  superior  to  triangular 
elements.  This  experiment  re  implements  the  primitive  model 
using  equilateral  elements  and  a  comparison  is  made  between 
the  results  of  both  implementations. 

A.  MCIZL  DESCRIPTION 

A  form  of  the  barotropic,  shallow  water,  primitive 
equations  developed  by  Phillips  (1959)  is  used  as  the 
governing  equation  set  for  this  model.  In  Cartesian 
coordinates  the  equation  set  is 
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The  finite  element  formulation  of  this  set  of  equations 
evaluates  the  height  and  velocity  components  at  the  same 
nodal  points.  This  is  an  important  consideration,  because 
the  other  models  (the  linearized  model,  see  Chapter  VI,  and 
the  vorticlty-divergence  model,  see  Chapter  III)  either 
stagger  the  dependent  variables  or  have  the  property  of  a 
staggered  formulation.  When  comparisons  are  made  between  the 
models,  it  is  this  lattice  structure  that  is  being  compared. 


This 

form 

of 

the  shallow  water 

equations  includes 

avity 

waves 

as 

a  solution.  Gravity 

waves  have  a  maximum 

phase  speed  of  about  300  meters/second.  When  the  correct 
time  step  is  calculated  using  equation  III-26,  a 
considerably  smaller  time  step  is  obtained  compared  to  the 
larger  time  step  permitted  in  the  vort icity-divergence 
formulation.  This  is  an  important  feature.  If  solutions  from 
all  models  are  equally  as  good,  the  best  formulation  would 
be  determined  using  the  computational  time  required  to 
produce  the  desired  forecast. 
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All  models  use  the  same  domain  structure.  In  fact,  the 
detrain  described  in  Chapter  III  was  patterned  after  the 
domain  implemented  hy  Kelley.  Again,  this  domain  simulates  a 
belt  around  the  earth,  with  cyclic  continuity  which 
eliminates  the  east-west  boundaries.  Rigid  boundary  walls 
exist  at  the  equator  and  at  30  degrees  north  latitude.  The 
domain  is  composed  of  a  12x12  point  mesh  and  subdivided  into 
the  right  triangular  elements  illustrated  in  Figure  9. 
Notice  that  the  grid  points  are  not  shifted  as  in  the 
equilateral  element  implementation  shown  in  Figure  6. 

The  following  boundary  conditions  are  imposed: 

1)  -  no  cross  channel  flow  at  the  latitude 

boundaries . 

2)  -  a  geostrophic  balance  at  the  channel  walls 

imposed  on  the  continuity  equation  V-l. 

This  model  has  a  simple  second  order  diffusive  term  in 
the  equations  of  motion  7-2  and  V-3.  However,  for  the 
purpose  of  evaluating  these-  different  f ormul  at- loirs ,  -this 
option  was  not  implemented  during  the  comparison  phase. 

Initial  conditions  consist  of  a  single  wave  in  the  x 
direction  and  a  half  wave  in  the  y  direction.  The  initial 
fields  for  the  three  dependent  variables  are  shown  in  Figure 
10.  The  maximum  zonal  wind  perturbation  of  5.5  meters/second 
is  superimposed  on  a  mean  zonal  flow  of  10  meters/second. 
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Figure  9  .  Domain  divided  Into  right  triangles. 
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Figure  10.  Initial  flelda  for  the  primitive  model.  Both  the  x  and  y 
axes  are  multiplied  by  10  Km. 
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This  cursory  description  of  the  primitive  model  mentions 
only  those  significant  features  that  will  weigh  heavily  in 
the  ccmparisions  later.  The  Galerkin  implementation  of  this 
model  is  similar  to  that  presented  in  Chapters  II  and  III 
and  the  system  of  equations  is  solved  using  a  Gauss-Seidel 
iterative  procedure.  Further  details  concerning  this 
primitive  model  are  given  in  Kelley  (1976). 

B.  P.ISDLTS 

This  experiment  involves  the  shape  of  the  elements.  As 
mentioned  previously,  Kelley's  implementation  subdivides  the 
domain  into  right  triangular  elements,  as  illustrated  in 
Figure  9.  Considerable  small-scale  noise  was  observed  by 
Kelley  in  the  48  hour  forecast  solution. 

The  transition  from  right  triangles  to  equilateral 
triangles  changes  the  size  of  the  domain.  With  right 
triangles,  the *x  and«y  grid  spacings  are  equal  (380  KM).  A 
12x12  grid  matrix  has  a  length  and  width  of  3600  KM.-  With 
equilateral  triangles,  the*x  and  ay  grid  spacings  are  no 
longer  equal.  Arbitrarily,  the  ay  grid  spacing  is  held 
constant  (300  KM)  and  a  new*x  grid  spacing  computed  by 

ax  *  Ay/ COS (30)  V-4 
A  12x12  grid  matrix  with  equilateral  elements  has  a  width  of 
3600  KM  and  a  length  of  4045  KM. 

m 

Figure  11  contains  the  48  hour  forecasts  produced  using 
both  types  of  elements.  The  three  dependent  variables  fields 
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Figure  11.  48  hour  forecast  comparison  from  the  primitive  model  using 

both  right  triangles  and  equilateral  triangles  zor  element  shapes. 
Arbitrary  perturbation  amplitude  of  5.5  m/s. 


are  compared  for  each.  The  small-scale  noise  that  Kelley 
observed  is  present.  The  three  fields  shew  varying  degrees 
cf  distortions,  which  are  especially  noticeable  along  the 
boundaries . 

Older  ( 1S81 )  found  that  th^oot  mean  square  error  (SKSE) 
was  reduced  22  percent  by  using  equilateral  shaped  elements. 
This  improvement  is  apparent  on  viewing  Figures  11b,  d  and 
f.  The  contours  are  smooth  and  the  boundaries  are 
noise  free.  Kelley  showed  excellent  treatment  of  wave 
propagation  by  this  primitive  model.  The  lowest  resolution 
grid  '6x6)  tested  by  Kelley  was  within  four  percent  of  the 
true  phase  velocity.  Changing  the  element  shape  had  no 
apparent  effect  on  the  phase  velocities. 

lecause  the  outcome  of  this  experiment  was  a 
significantly  improved  forecast  solution,  future  comparisons 
with  the  primitive  model  will  be  made  using  equilateral 
elements. 
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"he  previous  chapter  r emonst rated  how  the  snare  of  the 
element  oar.  significantly  improve  the  solution,  villiaxs  and 
7i ?n> i ewi cz  '  1951 }  used  differently  shaped  basis  functions 
on  a  linearized  eqraticr.  set  to  produce  excellent  solutions 
when  applied  to  the  gecstrophic  adjustment  protlerr. 

Spatial  staggering  of  dependent  variables  in  finite 
difference  formulations  has  given  rrucn  better  solutions  tc 
the  geostrophic  adjustment  process,  and  these  forms  are 
widely  used  in  meteorology.  Schoenstadt  '1930'  found  similar 
results  with  finite  element  formulations  with  piecewise 
lines-- basis  functions.  However,  staggering  nodal  points  is 
not  a  convenient  rethod  to  implement,  especially  in 
two-dimensions  with  irregular  boundaries,  so  alternative 
schemes  are  needed. 

The  implementation  of  the  alternative  scheme  introduced 
by  villiars  and  Zienkievicz  (1931)  are  presented  in  this 
chanter.  As  mentioned  above,  this  formulation  uses  different 
basis  functions  for  the  height  and  the  velocity  fields.  One 
o*’  the  basis  functions  is  piecewise  linear,  while  the  other 
Is  nie-ewise  constant,  as  is  illustrated  in  Figure  12  for  a 
one  dimensional  domain. 


f 


b)  Piecewise  constant  basis  function. 


Figure  12  .  Different  shaped  basis  functions. 
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A.  ICUATION  REFORMULATION 

The  primitive  form  of  the  shallow  water  equations 
presented  in  Chapter  V  (V-l.V-2  and  V-3)  is  used  to  derive 
the  linearized  equations  needed  for  this  experiment .  The 
velocity  equations  V- 2  and  V-3  remain  unchanged  and  a  linear 
basis  function  (7  .  )  is  used  to  approximate  the  u  and  v 

J 

variables. 

The  continuity  equation  V-l  is  linearized  as  follows: 


negligible 

where  f  is  the  average  geopotential  over  the  domain.  A 
piecewise  constant  basis  function  (V j  )  is  used  to 
approximate  the  geopotential.  This  linearization  is 
reasonable  in  this  case  because  the  Rossby  radius  of 
deformation  J'^/f.is  much  larger  than  ax  [see  Williams  and 
Zlenkiewicz  (1981)].  The  Galerkin  method  is  applied,  to.  thi s 
linearized  equation  set  using  a  piecewise  linear  test 
function  for  7-2  and  V-3,  and  a  piecewise  constant  test 
function  for  VI-1. 

The  piecewise  constant  basis  function  has  the  property 
of  displacing  the  geopotential  to  the  centroid  location  of 
the  elements,  which  should  give  the  same  effect  as 
staggering  the  grid  points,  as  seen  in  Figure  13.  The 
density  of  geopotential  data  in  the  domain  is  now  greater 
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Figure  13  .  Staggering  of  the  geopotential 
about  the  velocity.  The  o  symbol  are  the 
actual  grid  point  locations  and  the  •  symbol 
the  centroid  location  of  each  element. 


than  the  density  of  the  velocity  data.  Instead  cf  a 
geooo tent ial  value  for  each  node,  there  is  one  averaged 
value  ever  each  elerrent. 

The  final  form  of  the  Galerkin  equation  for  VI-1,  V-2 
and  V-3  after  performing  the  tire  differencing  is: 


>r<W  *  --2 

“TS-v  ■  -  •tU3<'w  *  “K'v'tac-v 

♦  VjUk'7j  *Vky  ,V1>  ~  f°Vj<7j  ,7i>]  7I"3 


n+1 . 


'PW  *  *  -  *tH3<ijyvi>  *  ’JjvS<'rj-'rk*-vi> 

'  f.“3<»yV1>]  VI-4 


This  linearized  set  of  equations  VI-2,  VI-3  and  VI-4  is 
solved  using  a  lauss  Seidel  iterative  procedure.  It  is  worth 
mentioning  that  the  coefficient  matrix  <¥.,¥.>  in  equation 
VI-2  has  all  non  zero  coefficients  equal  to  one,  since  the 
integration  of  the  inner  product  involves  piecewise 
constant  functions. 

This  equation  set  is  implemented  rather  than  the 
equation  set  V-l,  V-2  and  V-3  using  all  cf  the  existing 
primitive  model  code.  The  major  modification  involved  the 
way  that  the  geopotential  was  referenced.  With  an  average 
geopotential  over  each  element  Instead  of  a  value  at  each 
node  for  a  12x12  mesh,  there  are  2S8  geopotential  points 
versus  the  Iff  velocity  (u,v)  points. 
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E.  PZSDITS 

The  results  from  this  linearized  model  are  compared  to 
those  from  the  primitive  model.  The  initial  field  for  this 
experiment.  Figure  14  has  a  maiimium  perturbation  zonal  wind 
that  is  one  fifth  of  the  value  used  in  Chapter  V,  Figure  10. 
The  mean  zonal  wind  remains  10  meters/second.  The  48  hour 
forecast  solutions  for  each  field  are  compared  in  Figure  15. 

This  alternative  formulation  shows  some  promise, 
although  there  are  some  minor  perturbations  in  these 
contours  compared  to  those  in  the  primitive  solution.  No 
explanation  is  offered  for  this  small-scale  noise,  although 
possibly  the  other  formulations  presented  b7  Williams  and 
Zier.kiewicz  (1981)  would  improve  the  solutions. 
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VII.  V ORT I C IT Y~II  £N Cl  MCTZL  ZXPIRJEMZM 

•As  in  the  previous  two  chapters,  a  comparison  between 
two  models  will  be  presented.  The  results  from  the 
vorticity-divergence  model  are  compared  to  the  results  from 
the  primitive  model.  To  fully  exploit  the  differences 
between  the  models  performances  and  indicate  the  strengths 
ard  weaknesses  of  each  formulation,  three  domain  geometries 
and  three  initial  conditions  are  used.  All  solutions  are  at 
43  hour,  except  for  one  case  which  was  extended  to  36  hours. 

From  these  two  finite  element  formulations  sore 
additional  insight  is  obtained  concerning  the  execution  time 
required  as  the  grid  resolution  changes.  Lastly,  a  brief 
discussion  on  the  sensitivity  of  the  computational  technique 
is  miven. 

A.  TIST  TOMAINS  AM  INITIAL  CONDITIONS 

The  three  domain  geometries  used  in  the  model  evaluation 
are  illustrated  in  Figure  16.  All  domains  consist  of  a  12x12 
elerent  mesh  with  equilateral  shaped  elements  (156  grid 
points'  and  cyclic  continuity  is  imposed  on  the  east  and 
west  boundaries.  The  domain  has  dimensions  of  4045  KM  along 
the  x  axis  and  35?3  KM  along  the  y  axis. 

The  regular  domain  (Figure  16a)  has  a  uniform 
distribution  of  grid  points,  with  a  minimum  grid  spacing 
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along  the  x  axis  of  337  XM.  This  is  the  most  congenial 
domain  and  oroduces  the  test  results. 

The  smooth  domain  (Figure  15b)  has  a  smoothly  varying 
distribution  of  grid  points.  This  technique  allows  a  smooth 
variation  of  resolution.  It  was  developed  by  Older  ( 1361). 
who  showed  how  it  significantly  reduced  noise  at  all 
frequencies  compared  with  other  variable  grid  domains. 

The  degree  of  resolution  variation  is  acccrrpl ished  by 
selecting  an  appropriate  value  for  the  ratio  of  maximum 
stretch  to  minimum  shrink  along  both  axes.  The  experiments 
presented  in  this  chapter  use  a  2.5  ratio  for  both 
directions.  This  produces  a  grid  point  concentration  in  the 
right  center  of  Figure  16b  and  coarse  resolution  at  the  to? 
and  bottom  left  o?  the  domain.  The  minimum  grid  spacing 
along  the  x  axis  is  1SS  XM  in  the  dense  grid  point  area. 

The  third  domain  was  the  least  hospitable  geometry  for 
both  models.  The  abrupt  domain  (Figure  16c)  has  a  dense  grid 
point  concentration  on  the  left  and  coarse  spacing  cn  the 
right  of  the  domain.  The  minimum  grid  spacing  along  the  x 
axis  is  168  XM  in  the  area  on  the  left.  The  grid  spacings 
are  uniform  except  for  the  abrupt  change  along  the  center. 

Although  these  three  domains  are  simple  in  structure, 
they  are  adequate  to  test  both  models'  performance.  To 
further  enhance  some  differentiating  characteristics  between 
the  two  formulations,  three  initial  conditions  are  used.  Two 
have  previously  been  described  in  Chapters  7  and  71.  Their 
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rain  distinguishing  feature  is  the  arbitrary  perturbation, 
amplitude  'A?A)  magnitude.  The  nearly  linear  case  has  an  APA 
=  1.1  m/s  ccxpared  tc  5.5  r/s  for  the  xore  nonlinear  case. 
All  initial  conditions  have  a  10  m/s  rean  flow  component. 

Vith  the  nearly  linear  case  both  models  behave  well.  The 
introduction  of  the  more  nonlinear  initial  disturbance 
illustrated  the  boundary  and  computational  technique 
sensitivity.  The  third  initial  condition  is  the  nearly 
linear  initial  field  with  the  wave  length  equal  to  half  the 
dorain  length.  This  has  the  effect  of  producing  two  waves 
with  the  domain  . 

2.  TEST  CASE  COMPARISONS 

The  comparisons  between  models  are  divided  according  to 
the  domain  geometry.  The  primitive  model  has  three  fields; 
geopotential,  o  and  v.  The  vorticity  —  divergence  model  has 
seven  fields:  geopotential,  stream  function,  velocity 
potential,  u,  v,  vorticity  and  divergence.  The  geopotential , 
u  and  v  will  be  the  only  fields  used  for  the  comparison. 

1 .  Regular  Case 

The  regular  domain  geometry  (see  Figure  16a)  is  used 
'or  this  first  comparison.  The  initial  conditions  have  an 
A? A  *  1.1  m/s,  as  shown  In  the  contour  plots  ir.  Figure  17. 

The  48  hour  solutions  for  both  models  are  presented 
in  Figure  18.  As  anticipated,  all  of  the  solutions  have 
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Figure  17.  Initial  fields  for  both  the  primitive  and  vorticity — 
divergence  models  using  the  regular  domain  and  perturbation  amplitude 
of  1.1  m/s. 

88 


smoothly  contoured  fields  with  no  noticeble  noise  and  their 


phase  velocities  are  comparable. 

Vi th  the  nearly  identical  solutions,  the 
distinguishing  feature  between  the  models  is  the 
computational  time.  The  computational  time  is  determined  by 
the  size  of  the  time  step  each  model  is  allowed  to  use  as 
indicated  on  eauation  III-25.  For  this  domain  the  grid 
suacins  is  33?  Z* .  The  maximum  phase  velocity  possible  is 
different  for  each  model.  The  primitive  model  allows  gravity 
waves  so  c  =  320  m/s  whereas  the  vorticity-divergence  model 
filters  out  these  high  frequency  waves.  A  c  =  10  m/s  is  used 
for  this  formulation. 

Table  1 


Vcdel 

ox 

c 

at 

#  cf  steps 

CPU  units 

(IV ) 

(m/s) 

(sec) 

(iterations) 

(  sec) 

PE 

337 

300 

45S 

3?e 

145. 5 

V-E 

33? 

10 

13,755 

13 

33.0 

omoarisons  of 

the 

computational 

times  for 

a  43  hour 

forecast  between  the  primitive  model  and  the 
vorticity-divergence  model. 

Table  1  compares  the  results  of  the  computational 
times  for  both  models.  The  vortlcity  divergence  model 
produced  as  accurate  results  over  the  regular  domain  using 
22?  of  the  computational  time  needed  by  the  primitive  model. 
In  fact  from  geostrophic  reasoning  the  vort icity-divergence 
model  should  produce  better  forecasts  when  small  scale 
features  are  present. 
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2 .  Smooth  Case 


The  smooth  domain  geometry  (Figure  let)  is  used  in 
this  second  comparison.  The  initial  fields  shown  in  Figure 
12  have  a  disturbance  amplitude  of  1.1  r/'s. 

The  48  hour  solution  comparison  is  presented  in 
Figure  20.  All  fields  again  have  smooth  contoured  plots. 
Close  inspection  of  the  two  u  fields,  Figures  20c  and  d, 
shows  that  the  orimitive  u  field  has  small  hicks  along  some 
contours  and  a  weak  tilt  near  the  central  boundary  nodes, 
although  this  may  be  a  function  of  the  plotting  routine. 
Notice  the  good  symmetry  for  the  vort ic ity-di vergence  u 
field. 

As  in  the  regular  case,  this  smooth  experiment 
produced  two  acceptable  solutions.  Again,  the  computati cnal 
time  is  the  differentiating  criteria. 

Table  2 


odel 

AX 

(Kt* ) 

(m/s  ) 

at 

(sec ) 

*  of  steps 
( i  teratior.s ) 

CPU  units 
( sec) 

PI 

129 

300 

271 

838 

304.9 

-D 

199 

10 

8124 

21 

42 . 5 

Comparison  of  the  computatio nal  times  for  a  43  hour  forecast 
between  the  primitive  model  and  the  vorticity-divergence 
model  using  the  smooth  domain. 

Table  2  comuares  the  computational  times  for  both  models. 
The  vortici ty-dlvergence  has  a  53.8?  saving  of  CPU  time. 
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Figure  19.  Initial  fields  for  both  the  primitive  and  vorticity- 
divergence  models  using  the  smooth  domain  and  perturbation  amplitude 
of  1.1  m/s. 


Ir.  an  effort  to  contrast  the  computational  accuracy 
cf  the  models,  this  experiment  is  repeated  using  the  "ore 
nr.iir.ear  initial  case,  Figure  21  shews  the  initial  fields 
with  A? A  =  5.5  rr/s.  This  larger  initial  disturbance  is 
reflected  ir.  the  greater  geopotential  amplitude  and 
magnitudes  of  the  contour  lines. 

Figure  22  shows  the  43  hour  solution  comparison.  As 
in  the  more  linear  case  presented  above,  all  o?  the  plots 
have  smooth  contours  with  no  noticeable  noise.  The 
vort ici ty-di vergence  geopotential  field.  Figure  22b,  has  the 
ridge  extending  farther  north,  and  flattening  of  the 
southernmost  contour,  than  does  the  primitive  geopotential 
field,  Figure  22a.  The  mean  geopotential  heights  for  both 
models  have  also  increased.  The  primitive  mean  geopotential 
is  new  49280  gum  and  the  vortici ty-divergence  geopotential 
is  49622  gpm. 

These  two  discrepancies  indicate  that  the  boundaries 
are  not  handled  accurately.  Turing  the  implementation  of  the 
vortioity-dlvergence  model,  treatment  of  the  boundaries  was 
the  most  troublesome  phase.  The  vortici ty-divergence 
formulation  is  a  complex  equation  set  and  time  limitations 
restricted  further  investigation  of  more  sophisticated 
boundary  conditions. 

This  same  initial  condition  is  now  extended  to  a  56 
hour  forecast,  which  is  shown  in  Figure  23.  The  mean 
vcrticity-divergence  geopotential  increased  tc  4S600  gpm. 
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Figure  21.  Initial  fields  for  both  the  primitive  and  the  vorticity- 
divergence  models  using  the  smooth  domain  and  perturbation  amplitude 
of  5.5  m/s. 
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Figure  22.  ^6  hour  forecast  comparison  between  the  primitive  model 

and  the  vorticity-dlvergence  model  using  the  smooth  domain.  (APA  -  5*5 
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Figure  23*  96  hour  forecast  comparison  between  the  primitive  model 

and  the  vorticity-dlvergence  model  using  the  smooth  domain.  (APA  - 


whereas  the  irean  primitive  geopotential  remained  constant. 
All  fields  have  smooth  contours.  Close  inspection  of  both  u 
fields,  Figure  22 c  and  i,  shows  a  slight  skewing  of  the 
contours  along  the  central  channel  grid  points.  The 
vorticity  divergence  u  field  has  the  most  pronounced 
deviations.  A  hypothesis  for  this  skewing,  explained  further 
below,  is  that  it  is  caused  by  the  corputat ional  technique 
employed.  The  relaxation  schemes  are  extremely  sensitive  and 
fine  tuning  of  the  relaxation  coefficient  would  have 
required  more  time  than  was  available. 

Table  3 


mod  el 

AX 

(m/s) 

iA  &* 
CD  <-#■ 
O 

#  of  steps 
( iterations ) 

CPU  units 
(  sec) 

py 

199 

300 

271 

1276 

5S4.0 

V-D 

139 

10 

8124 

42 

31.0 

Comparison  of  the  computational  tires  for  a  S€  hour  forecast 
between  the  primitive  model  and  the  vortici ty-divergence 
model  using  the  smooth  domain. 


Table  3  shows  the  comparison  between  both  models  for 
the  56  hour  forecast.  A  savings  of  55  percent  is  realized 
with  the  vorticity  divergence  model. 

The  above  experiment  points  out  the  two  areas  where 
the  vorticity  —  divergence  model  is  presently  weak,  the 
increase  of  the  geopotential  and  the  sensitivity  of  the 
relaxation  coefficients.  3oth  these  weaknesses  can  be 
improved  and  are  not  a  result  of  the  formulation  but  of  the 
implementation.  At  oresent,  their  Influence  is  not  detected 
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er-ept  in  extremely  Ion,?  forecasts,  such  as  in  the  96  hour 
example.  Further  experiments  ray  still  he  required  if  they 
demonstrate  significant  differences  in  the  solutions  c:'  the 
two  models. 

"he  last  experiment  on  the  smooth  dorain  uses  the 
mere  linear  initial  condition  A?A  =  1.1  m/s,  tut  its  wave 
length  is  divided  by  two,  so  that  two  shorter  waves  are 
propagated  through  the  channel.  The  initial  conditions  are 
shown  in  Figure  24.  Decreasing  the  wave  length  has  the  same 
effect  as  decreasing  the  density  of  grid  points.  In  this 
case  six  grid  points  are  used  to  describe  the  wave  structure 
versus  the  12  used  in  the  previous  cases. 

The  48  hour  comparisons  are  shown  in  Figure  25.  As 
in  all  previous  cases,  a  computational  time  saving  of  54*  is 
gained  with  the  vorticity-divergence  model,  With  fewer  grid 
points  describing  the  wave  structure,  more  small-scale  r.oise 
is  introduced  into  the  solution.  Comparing  the  primitive 
geopotential  field.  Figure  25a,  to  the  initial  geopotential, 
Figure  24a,  shows  a  dampening  of  the  wave  amplitude,  whereas 
the  vorticity-divergence  geopotential.  Figure  25b, 
correlated  well  with  the  initial  geopotential. 

The  high  frequency  noise  is  evident  on  both  u 
fields.  Figures  25c  and  d.  The  primitive  model  u  field  is 
poorly  defined  along  the  boundary  and  becomes  irregular  over 
the  interior  grid  points. 
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Figure  25.  hour  forecast  comparison  between  the  primitive  model 
and  the  vorticity-divergence  model  using  the  smooth  domain.  Two  waves 
are  embeded  in  the  flow  and  APA  ■  1.1  m/s. 
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3.  Abrupt  Case 

The  abrupt  case  comparison,  uses  the  abrupt  domain 
geometry  (see  Figure  16c).  This  grid  point  configuration  is 
used  to  *urther  illustrate  the  effects  of  spatial 
resolution.  The  previous  case  using  the  smooth  domain  and 
half  wave  length  introduced  noise  into  the  solution,  but  the 
spatial  resolution  changed  slowly  and  gradually. 

Consider  the  transition  necessary  in  an  operational 
model,  where  the'luiury  of  iiavlnp  a  long  smooth  transition 
Into  the  region  of  high  resolution  may  not  be  possible.  The 
abrupt  domain  is  an  example  of  the  results  obtained  when 
spatial  resolution  is  decreased  rapidly. 

The  initial  fields  are  shown  in  Figure  26.  The  more 
linear  case,  A?A  *  1.1  m/s,  is  used  to  eliminate  effects  due 
to  the  Initial  field,  so  that  only  the  effects  due  to  the 
grid  geometry  are  seen. 

The  48  hour  comparisons  are  shown  in  Figure  27.  3oth 
solutions  are  affected  by  this  geometry,  but  the  primitive 
solutions  are  totally  disorganized  and  unacceptable. 

Table  4  shows  the  comparison  of  computational  times 
for  both  models  for  a  48  hour  forecast.  An  353;  savings  in 
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A!  INITIAL  0  FIELD  ICPM) 


8)  INITIAL  U  FIELD  (M/S) 


3 


Figure  26.  Initial  fields  for  both  the  primitive  and  vorticity- 
divergence  models  using  the  abrupt  domain  and  perturbation  amplitude 
of  1.1  m/s.  Note  that  the  element  shapes  are  not  equilateral  triangles 
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Figure  27.  48  hour  xot*>  ,as«.  :omparison  between  the  primitive  model 
and  the  vorticity-divergence  model  using  the  abrupt  domain.  (APA  ■ 
1.1  m/ s ) . 
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CPU  time  is  obtained  using  the  vorticity-divergence  model. 
Met  only  is  there  a  computational  savings  with  the 
vorticity-divergence  model,  but  the  solution  is 
significantly  better  than  the  primitive  model. 


Table  4 


Model 

A  X 

'ry ) 

(m/s ) 

At 

(sec ) 

w  of  steps 
(iterations ) 

CPU  units 
(sec) 

PF 

165 

300 

228 

756 

368.9 

V-E 

166 

10 

6858 

25 

55.7 

Comparison  of  the  computaional  times  for  a  48  hour  forecast 
between  the  primitive  model  and  the  vcrticity-divergence 
model  using  the  abrupt  domain. 

C.  COMPUTATIONAL  SLNSITIVITT 

This  section  will  offer  an  explanation  for  the  skewing 
0^  the  contours  in  the  56  hour  forecast  solut ion- over  the  - 
smooth  domain  -Figure  23).  As  mentioned  oreviousiy,  there 
was  rot  enough  time  available  for  fine  tuning  the 
overrelaxation  coefficient,  which  is  used  while  solving  the 
system  of  equations.  The  overrelaxation  coefficient  may  be 
very  sensitive  and  small  changes  can,  on  occasion,  radically 
change  the  rate  of  convergence.  The  optimal  value  of  the 
overrelaxation  coefficient  depends  on  the  specific  form  of 
the  coefficients  of  the  equation  and  the  error  distribution. 

The  equation  set  to  be  solved  consists  of  three 
equations  and  each  equation  required  its  own  relaxation 
coefficient.  When  solving  tb  equations  over  the  regular 
domain,  the  entire  system  is  well  behaved  and  an  optimum 
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relaxation  coefficient  is  easily  determined.  However,  as  the 
dorain  <?ecmetry  changes,  the  system  of  equations  ic  net 
converse  as  raoidly  and  the  relaxation  coefficient  requires 
*urther  fine  tuning. 

The  9€  hour  forecast  uses  the  smooth  domain.  The 
mid-latitudinal  grid  points  are  compacted,  creates  a  denser 
belt  in  the  middle  of  the  channel.  The  coefficients 
originally  computed  using  the  regular  domain  need 
adjustments  to  properly  solve  the  equations. 

To  illustrate  the  significance  for  fine  tuning  the 
relaxation  coefficient,  consider  the  series  cf  plots  in 
Figure  29.  The  vorticity  divergence  equation  set  can  he 
simplified  by  assuming  the  flow  is  non-divergent ,  so  that 
only  the  vorticity  equation  needs  to  be  solved.  Figure  '29a 
is  the  46  hour  vorticity  field  using  this  equation  over  the 
regular  domain  with  an  overrelaxation  coefficient  of  1.3. 
The  field  is  well  defined  with  smooth  contours. 

Figure  28b  is  similar  to  the  case  in  Figure  28a  except 
that  the  smooth  domain  is  used.  Notice  the  V-shaped  kink  in 
the  uattern  with  a  steeper  slope  in  the  upper  half.  This 
increased  bias  in  the  upper  half  is  caused  by  relaxing  the 
field  in  the  same  direction  during  each  pass  over  the 
domain.  V/hen  the  direction  is  reversed  after  each  pass,  the 
exaggerated  bias  in  the  upper  half  disappears,  as  is  seen  in 
Figure  28c.  However,  the  7-shaped  kink  is  not  eliminated. 
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Figure  28.  Computational  sensitivity  using  the  48  hour  vorticity  fields. 
Fig.  28  A)  the  48  hour  forecast  using  the  regular  domain,  28  B)  the  same 
forecast  using  the  smooth  domain  relaxing  in  one  direction,  28  C)  same  as 
28  B)  except  alternate  the  direction  of  relaxation  after  each  pass  over 
the  domain,  and  28  D)  same  as  28  C)  except  the  relaxation  coefficient  was 
fine  tuned  form  1.3  to  1.297. 
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Varying  the  relaxation  coefficient  frcrr  1.3  to  1.2S? 
produces  the  much  improved  solution  in  Figure  28d.  If  the 
relaxation  coefficients  for  the  other  two  equations  could 
also  be  fine  tuned,  it  appears  that  improved  solutions  would 
result.  There  are  also  other  relaxation  techniques  available 
that  have  potential  for  improving  the  solution  while  also 
converging  at  a  faster  rate.  Some  of  these  techniques  will 
be  tested  in  the  futuer  using  this  model. 
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VIII.  CONCLUSIONS 


This  research  investigated  different  finite  element 
f orrulations  for  the  shallow  water  equations.  The 
two—  dimensional  domain  was  a  channel  which  simulated  a  belt 
around  the  earth.  Analytic  initial  conditions  were  used  to 
simplify  the  comparisons.  Two  formulations  were  examined; 
one  using  different  shaped  basis  functions  and  the  other 
using  a  different  form  of  the  equation.  lach  was  compared  to 
the  urirritive  form  of  the  shallow  water  equations  that  was 
developed  by  Kelley  (19?S). 

The  use  of  equilateral  shaped  elements  which  was 
suggested  by  T>.  M.J.F.  Cullen  significantly  improved  the 
solutions  compared  to  Kelley's  model,  which  originally  used 
right  triangles  as  basis  functions.  Most  of  the  other 
studies  in  this  thesis  used  the  equilateral  triangles. 

Williams  and  Zlenhiewicz  (1981)  suggested  the  use  of 
piecewise  linear  basis  functions  for  the  velocity  field  and 
piecewise  constant  functions  for  the  height  field.  This 
formulation  was  tested  with  a  linearized  continuity 
equation.  The  results  were  poorer  than  those  obtained  with 
Kelley's  model. 

Most  of  the  effort  in  this  thesis  was  devoted  to 
lmplementating  and  testing  a  vorticity-di vergeace  model 
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similar  to  the  ones  developed  by  Staniforth  and  Mitchell 
fl977)  and  Cullen  and  Hall  (1S?9).  Several  tests  were 
presented  which  cotrpare  this  formulation  with  Kelley's 
model.  It  was  found  that  this  model  executes  approximately 
one  order  of  magnitude  faster  than  does  the  primitive 
formulation  used  by  Kelley.  Secondly,  as  the  spatial 
resolution  between  grid  points  decreases,  this  formulation 
oroduces  a  solution  that  is  far  superior  to  the  primitive 
form.  A  disadvantage  is  its  computational  sensitivity,  which 
requires  fine  tuning  in  solving  the  elliptic  equations  for 
certain  geometries.  It  also  requires  25  percent  more 
computer  storage,  due  to  the  more  complex  equation  set  and 
the  additional  variables  that  are  treated. 

Implementation  of  finite  element  models  is  not  easy. 
Fowever,  there  is  a  lot  of  generality  and  redundancy 
irbedded  in  the  computations.  Versatile  modules  were  written 
which  significantly  reduced  the  effort  in  implementing  the 
vorticity-divergeace  model. 

Further  research  is  suggested  using  this  finite  element 
formulation.  It  has  accurate  phase  propagation,  is  able  to 
handle  variable  grid  geometry,  reduce  the  small-scale  noise 
and  decrease  the  model's  execution  time.  Specifically  more 
advanced  methods  of  solving  the  elliptic  equations  should  be 
investigated,  finally,  the  formulation  should  be  tested  with 
small-scale  forcing,  where  its  advantages  should  be  most 
evident . 
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